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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC)  UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  (met¬ 
ric)  units  as  follows: 

Multiply  By  To  Obtain 

feet  0.3048  metres 

inches  2.54  centimetres 
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WAVERTF:  WAVE  PROPAGATION  OVER  OBSTACLES 


AND  IRREGULAR  TOPOGRAPHY 
A  USER'S  MANUAL 

PART  I:  MATHEMATICAL  FORMULATION 

1.  In  this  study,  a  hybrid  element  method  is  developed  for  solving  the 
progagation  of  linear  water  waves  over  a  finite  near  field  involving  irregular 
obstacles  and  bathymetry.  In  the  far  fields  water  depths  are  assumed  to  be 
constant  but  not  necessarily  equal.  A  geometrical  configuration  of  the  prob¬ 
lem  is  illustrated  in  Figure  1 .  The  problem  has  been  studied  by  numerous  in¬ 
vestigators  over  the  years  (Bai  and  Yeung  1974;  Black,  Mei,  and  Bray  1971; 

Chen  and  Mei  1974;  Dean  and  Ursell  1959;  Lamb  1945;  Lee,  Ayer,  and  Chiang 
1980;  Miles  1967;  Mynett,  Serman,  and  Mei  1979;  and  Newman  1965a, b).  How¬ 
ever,  their  problems  are  often  oversimplified,  and  the  solutions  are  not  ef¬ 
fective  in  application.  Effective  application  of  these  problems  can  be  found 
in  Chen  (1984)  where  the  solutions  are  similar  to  those  in  Part  I  of  this 
report . 
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Boundary  Value  Problem 


2.  Let  (x,z)*  be  Cartesian  coordinates  with  z  =  0  representing  the 
undisturbed  water  surface  and  with  the  upward  direction  being  the  positive 
z-axis,  as  shown  in  Figure  1.  The  fluid  motion  in  a  linear  wave  field  is  as¬ 
sumed  to  be  two-dimensional  such  that  the  conservations  of  mass  and  momentum 
are  written  as 


J 


(2) 


where 

Uj  =  J-component  flow  velocity 

J=1,2  =  x  and  z  components,  respectively 
t  =  temporal  coordinate 
p  =  pressure 
p  =  water  density 
g  =  gravitational  acceleration 

3.  In  Equation  2  the  friction  term  is  analogous  to  Heaps  (1969)  and  is 
modeled  by  -eUj  which  is  linearly  proportional  to  the  flow  velocity.  The 
coefficient  e  is  generally  a  spatial  complex  function,  indicating  spatial 
variability  and  phase  difference  from  Uj  .  A  more  detailed  explanation  of 

e  will  be  given  later. 

4.  For  convenience,  a  potential  function  $  is  defined  such  that 


i£  .  2 

dt  ~  p 


*  gz 


(3) 


Since  the  wave  is  assumed  to  be  sinusoidal  in  time  with  radian  frequency  w  , 
we  may  separate  temporal  and  spatial  dependencies  by 


*  For  convenience,  symbols  and  abbreviations  are  listed  in  the  Notation 
(Appendix  C). 
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Uj(x,z,t) 

4>(x,z,t) 

t(x,z,t) 


Uj(x,z)' 

$(x, z ) 
n(x,z) 


exp  (-iwt) 


where 


C  =  the  free  surface  displacement 
U  =  spatial  part  of  flow  velocity 
•  =  spatial  part  of  the  velocity  potential  function 
i  =  VT 

n  =  free  surface  displacement 
Substituting  Equations  3  and  4  into  Equation  2,  we  have 


u.  =  X  !*- 
j 


where 


1  +  — 
to 

Substituting  Equations  4  and  5  into  Equation  1,  we  have 

V-XV*  =  0  ( 

where  v  is  the  two-dimensional  gradient  operator. 

5.  At  the  free  surface  z  =  c  ,  the  atmospheric  pressure  is  taken  as 
the  reference  pressure  p  =  0  ;  therefore,  Equations  3  and  4  give  the  free 
surface  displacement  in  terms  of  *  . 


c  =  -  s  it  or  ■» s  r  * 


at  z  =  0 


The  linearized  free  surface  boundary  condition  is  then  obtained  from  the  kine¬ 
matic  condition  and  Equation  8  to  become 


-  r  * 2  0 


at  z  =  0 


6.  At  the  solid  boundaries  along  the  bottom  and  the  stationary  obsta¬ 
cles,  an  absorbing  boundary  condition  is  used. 
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at  z  =  -h  and  3B 


(10) 


where  n  is  the  unit  normal  vector  outward  from  the  water  region,  o  can  be 
interpreted  as  the  absorption  coefficient.  If  the  normal  velocity  vanishes, 
then  a  =  0  . 

7.  At  x  +  ±®  in  the  far  fields,  the  radiation  conditions  are  re¬ 
quired  for  the  reflected  and  transmitted  waves  and  ,  respectively,  to 
ensure  a  unique  solution. 

+ 

a  ♦  + 

—  ±  ik“«"  =0  at  x  ♦  i®  (11) 

where  k*  and  k~  are  the  wave  numbers  of  the  propagating  modes  in  the  far 
fields.  The  superscripts  +  and  -  from  here  on  are  referred  to  as  the 
quantity  in  the  region  A+  and  A"  (see  Figure  1). 

8.  Therefore,  a  boundary  value  problem  is  established  with  Equation  7 
as  the  governing  equation  and  Equations  9,  10,  and  11  as  the  boundary  condi¬ 
tions.  Note  that  if  there  is  no  friction,  e  =  0  and  \  =  1  .  Equation  7 
reduces  to  a  Laplacian  and  Equations  9,  10,  and  11  to  their  counterparts  in 
the  usual  formulation  of  a  linear  wave  problem. 


Variational  Principle  and  Hybrid  Element  Approximation 


9.  A  hybrid  element  method  is  employed  to  oolve  the  boundary  value 
problem.  The  water  domain  is  divided  into  three  regions,  A  ,  A+  ,  A"  ,  as 
shown  in  Figure  1 .  A  conventional  finite  element  approximation  with  a  nodal- 
type  element  is  used  in  the  near  field  A  ,  and  analytical  solutions  with 
unknown  coefficients  are  used  to  describe  the  far  fields  A"  and  A+  .  A 
variational  principle  with  a  proper  functional  is  established  such  that  the 
matching  conditions  are  satisfied  at  the  common  boundaries  of  the  near  and  far 
fields.  Calculations  are  thus  localized  in  the  near  field  for  the  nodal  and 
coefficient  unknowns. 

10.  The  variational  principle  for  the  boundary  value  problem  requires 
that  the  following  functional  n  be  stationary  with  respect  to  the  arbitrary 
first  variation  of  ♦  and  •±  . 
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f f  i/2x(y«)2 

A 


dA 


/  Ti  “2  ds  -  /  10,2  ds 
3F  3B 


In  Equation  12,  «o(x,z)  is  the  spatial  part  of  the  incident  wave  and  is 

written  as 


« 

o 


(x,z) 


iga  cosh  k*(2  +  h+) 

^ - - exp  ('ikox) 

cosh  k  h 
o 


(13) 


where 


aQ  =  incident  wave  amplitude 
kQ  =  wave  number  of  the  propagating  mode 
Naturally,  the  incident  wave  is  coming  from  x  =  +®  and  propagating  to  the 
left.  Construction  of  the  functional,  Equation  12,  and  proof  of  the  equiva¬ 
lency  between  the  variational  principle  and  the  boundary  value  problem  are 
achieved  by  procedures  similar  to  those  given  by  Chen  and  Mei  (1974)  and  are 
not  given  in  this  report. 

11.  Triangular  elements  with  linear  shape  function  are  employed  to  sub¬ 
divide  the  near  field.  In  the  far  fields,  the  friction  terms  are  usually  of 
minor  concern  in  practice  and  are  omitted.  The  bottoms  are  also  assumed  to  be 
constant  depths  h*  and  impermeable  such  that  a  =  0  .  Therefore,  the 
analytical  solutions  of  the  outgoing  waves  in  the  far  fields  can  be  expressed 
as  follows: 


**  =  - 


iga. 


h 


.  cosh  k  (z  ♦  h”) 

R  - t— r -  exp  (±ik~x) 

ffl  cosh  k±h±  m 

m 


(14) 


where  the  coefficients  fr  are  unknown  constants  to  be  determined  and  are  the 


reflection  and  transmission  coefficients  when  m  =  0  .  k  and  k~  (m  >  1) 

o  m 

are  the  wave  numbers  of  the  propagating  and  evanescent  modes,  respectively, 

p 

and  are  determined  from  the  dispersion  relation  w  =  gk  tanh  kh  .  In  the 
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computation,  the  number  of  the  term  m  in  Equation  14  is  usually  truncated 

+  —4 

when  the  value  of  exp  (±ik"x)  is  the  order  of  10 

m 

12.  Using  the  conventional  finite  element  approximation  in  the  near 
field  and  Equation  14  in  the  far  fields,  one  can  deduce  a  set  of  simultaneous 
linear  equations  for  the  nodal  and  coefficient  unknowns  upon  extremization  of 
the  functional.  The  stiffness  matrix  is  symmetric  and  stored  in  a  packed  ma¬ 
trix  form  of  semibandwidth.  Gaussian  elimination  is  employed  for  a 
solution.  Computations  are  very  efficient. 

Free  Surface  and  Wave  Forces 


13.  Once  the  solution  is  obtained  for  the  velocity  potential  •  in  the 
near  field,  and  the  coefficient  constants  R*  in  the  far  fields,  the  free 
surface,  c  or  n  ,  is  given  by  Equation  8.  The  hydrodynamic  pressure  pd 
is  given  by  Equation  3  or 


Pd  =  "  p  ft  =  pi“*  exp  (“iwt) 


In  Equation  15  the  hydrostatic  pressure  is  excluded.  The  friction  is 


given  in  Equation  2  by  -eUj  or 


Tj  =  -  peUj  =  -pel  exp  (-iu>t)  (16) 

The  wave  forces  and  moment  on  the  obstacle  are  basically  the  integration  of 
the  hydrodynamic  pressure  and  friction  stresses  along  the  boundary  of  the  ob¬ 
stacle.  Therefore,  along  the  solid  boundary,  the  total  forces  Fj  and  the 
moment  Mj  in  the  positive  direction  perpendicular  to  the  (x,z)  plane  con¬ 
tributed  from  hydrodynamic  pressure  and  friction  are 


Fj  =/  <Pdnj  -  tj)  OS 

as 

fiy  =  f  f  x(pdn  -  x)  dS 


where 


n  =  (n  ,  n_)  =  unit  normal  outward  from  the  solid  boundary  3S 

X  z 
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\ 


r  =  distance  vector  from  the  reference  center  for  calcu¬ 
lating  moment 

T  =  (T. 


Tz> 


as  given  in  Equation  16. 

14.  From  Equations  8  and  15  through  18,  the  free  surface  n  ,  the  hy¬ 


drodynamic  pressure  , 
moment  My  are  thus  calculated 


the  friction 


TJ 


the  total  forces  F 


j 


and  the 


Friction  Terms 


15.  In  light  of  the  classical  solution  by  Stokes  for  an  oscillating 
flat  plate  in  laminar  flow  (Schlichting  1968),  in  which  the  retarding  force  is 
proportional  to  the  maximum  velocity  of  the  plate  and  has  a  phase  difference 
of  ir/4  from  the  velocity  of  the  plate,  e  is  expressed  as 

e  =  6  VvST  exp  (it)  ( 19) 

where  the  friction  coefficient  B  and  the  phase  difference  y  are  presumed 
to  be  dimensionless  constants  which  may  vary  spatially.  The  inclusion  of  a 
and  y  in  Equation  19  as  the  empirical  constants  is  intended  for  simplicity 
to  account  for  the  real  situation  in  most  engineering  applications  where  the 
boundary  and  bottom  are  not  smooth  and  flat.  For  the  Stokes  solution  men¬ 
tioned  above,  B  =  1  and  y  =  -x/4  .  Substituting  Equation  19  into 
Equation  6 


X  =  - - -  (20) 

'  *  uy/i  «‘T 

Since  Equation  19,  strictly  speaking,  is  not  a  universal  expression  for  the 
friction  coefficient,  some  other  expression  for  c  might  also  work  well.  The 
key  point  for  assessing  the  expression  of  e  will  be  simplicity,  efficiency, 
and  accuracy  in  application. 

Results  for  Idealized  Problems 


16.  Numerical  results  are  presented  for  the  propagation  of  water  waves 
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over  three  different  geometrical  configurations.  Very  little  is  known  about 
a  in  Equation  10,  and,  as  usual,  a  =  0  is  taken  in  all  of  the  following 
calculations. 

17.  Case  I  is  a  finite  step,  and  its  finite  element  network  is  shown  in 
Figure  2.  The  calculated  result  for  6  =  0  ,  as  shown  in  Figure  3,  agrees 


well  with  the  long  wave  theory  (Lamb  1945)  in  the  long  wave  range  and  the 
result  of  Newman  (1965a)  in  the  short  wave  range.  Note  that  Newman's  result 
is  from  the  case  of  an  infinite  step  where  h+  =  «  . 

18.  Case  II  is  a  fixed,  semi-immersed  circular  cylinder  in  a  finite 
water  depth.  Its  finite  element  network  is  shown  in  Figure  4.  The  calcu¬ 
lated  results  are  shown  in  Figures  5#  and  6  and  agree  fairly  well  with  the 
experimental  data  and  the  theoretical  results  from  Dean  and  Ursell  (1959), 


*  A  table  of  factors  for  converting  non-SI  units  of  measurement  to  SI  (metric 
units  is  presented  on  page  3. 
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Figure  4.  Finite  element  network  of  a  fixed,  semi-immersed 
circular  cylinder  in  a  finite  water  depth 


particularly  the  results  for  0  =  10  .  In  Figure  6,  the  horizontal  and 

vertical  force  coefficients  fx  and  fz  are  defined  as  |F  j  /(pga  a)  and 

\  x/max  ° 

(F  \  /(pgaQa)  where  ( F  I  and  (Fz)  are  the  maximum  forces  in  the 

z/max  \  x/max  \  z/max 


horizontal  and  vertical  directions,  respectively.  In  the  computer  program, 
the  forces  and  moment  coefficients  are  actually  defined  as 
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Figure  5.  Comparison  of  reflection  and  transmission  coefficients  R+  and  Rj 

the  fixed,  semi -immersed  circular  cylinder 


Figure  6.  Comparison  of  horizontal  and  vertical  force  coefficients  fx  and 

the  fixed,  semi -immersed  circular  cylinder 


f  = 


m.  = 


(Fx)max 

PgaoL 

(Fz) 

\ _ /max 


Pga0L 


Wmax 


pga0L 


(21) 


(22) 


(23) 


where 


y  •'/max 


is  the  maximum  moment  of 


in  Equation  18,  and  L  is  the 


length  scale  for  the  force  normalization  which  is  the  SCL(I)  in  the  computer 
program  given  in  Appendix  B.  L  =  a  is  used  in  this  calculation. 

19.  Note  that  the  theoretical  results  of  Dean  and  Ursell  (1959)  from 
the  assumption  of  potential  flow  in  an  infinite  water  depth.  In  the  calcula¬ 
tions,  the  friction  is  applied  only  on  the  solid  boundary  and  the  bottom  in 
the  near  field. 

20.  Case  III  is  a  rectangular  trench.  The  calculated  results  for  $ 

=  0  ,  as  shown  in  Figure  7,  agree  quite  well  with  the  theoretical  and  calcu¬ 
lated  results  and  the  experimental  data  by  Lee,  Ayer,  and  Chiang  (1980). 


LEGEND 


OO  O  O  PRESENT  RESULTS  USING  HYBRID  ELEMENT  METhO 


Figure  7.  Comparison  of  transmission  coefficient 
R~  of  a  rectangular  trench 
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PART  II: 


INPUT  DATA  PREPARATION  AND  STRUCTURE 


General 


21.  A  computer  code  WAVERTF  program  (listed  in  Appendix  B),  based  on 
the  hybrid  element  method  as  briefly  described  in  Part  I,  is  described  in  de¬ 
tail  here.  The  computer  code  is  implemented  using  FORTRAN  V  or  FORTRAN  Ex¬ 
tended  IV  on  the  CYBER  System.  To  use  standard  FORTRAN  IV,  changes  to  the 
program  must  be  made.  The  PARAMETER  statements  must  be  taken  out  and  the  val¬ 
ues  introduced  into  the  program  separately.  These  values  must  also  be  placed 
in  the  DIMENSION  statements.  Metric  (SI)  units  are  used  for  all  the  physical 
quantities  in  the  program.  If  other  units  are  used,  conversion  factors  must 
be  implemented,  mainly  the  viscosity  (VISCO)  and  the  gravitational  constant 
(G)  on  card  00130.  Note  that  a  unit  for  mass  is  immaterial  in  the 
calculation. 

Finite  Element  Grid 


22.  Some  guidelines  for  correctly  making  finite  element  grids  to  run 
the  computer  program  WAVERTF  are  given  below.  For  clarity,  we  shall  use  the 
case  of  a  fixed  semi-immersed  circular  cylinder  in  a  finite  water  depth,  shown 
in  Figure  4,  as  an  example. 

a.  The  incident  wave  is  from  x  =  +®  .  Therefore,  the  far  field 
at  the  right  hand  side  A+  will  be  the  reflected  region,  and 
at  the  left  hand  side  A“  will  be  the  transmitted  region. 
Between  the  two  far  fields  will  be  the  near  field  A  .  For 
convenience,  the  origin  of  the  Cartesian  coordinates  should  be 
chosen  somewhere  at  the  free  surface  in  region  A. 

b.  Before  starting  the  subdivision  of  the  near  field  domain  A 
into  elements,  the  outer  boundaries  of  the  A  domain,  3A+ 
and  a A-  as  shown  in  Figure  4,  must  be  established  to  accom¬ 
modate  the  obstacles  and  irregular  bottom.  The  outer  bound¬ 
aries  must  be  parallel  to  the  vertical  direction  (z-axis). 

c.  Triangular  elements  are  then  used  to  subdivide  the  domain  A  . 
For  best  results  the  grid  size,  although  arbitrary,  should 
change  gradually;  and  the  horizontal  extent  of  each  grid  size 
is  generally  not  greater  than  one-tenth  of  the  incident 
wavelength. 

d.  Numbering  of  the  nodes  requires  several  procedures. 

(1)  For  the  nodes  along  8A”  ,  the  numbering  should  start 
from  the  node  at  the  bottom  as  1  and  sequentially 
increase  along  3A~  to  the  node  at  the  free  surface. 
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(2)  For  the  nodes  along  3A+  ,  the  numbering  should  be  such 
that  the  node  at  the  free  surface  is  equal  to  the  total 
number  of  the  nodes  and  sequentially  decreases  along 

3A+  to  the  node  at  the  bottom. 

+ 

(3)  For  the  nodes  not  at  3A~  ,  the  numbering,  although  ar¬ 
bitrary,  should  be  such  that  the  maximum  value  of  the 
values  of  the  nodal  number  difference  in  each  element 
is  as  small  as  possible. 


Input  Data  Structure 

23.  To  run  the  computer  program  WAVERTF  (which  is  listed  along  with  the 
card  numbers  in  Appendix  B)  the  user  must  supply  input  data  for  a  specific 
problem.  The  detailed  requirements  for  these  user  inputs  are  given  in  this 
section.  An  example  problem  is  described  in  Part  III  along  with  a  listing  of 
the  corresponding  input  data. 

24.  There  are  only  two  PARAMETER  statements,  two  DATA  statements,  and  the 
subroutine  DATAIN  in  the  program  that  need  to  be  modified  for  each  new  problem. 

PARAMETER  and  DATA  Statements 


25.  The  numerical  values  in  the  following  four  statements  in  the  program 
need  be  modified  for  each  new  problem: 

00030  PARAMETER(NELE= 1 16 ,NN0D=8l ,NBD= 1 3,NKR=6 ,NKT=6 ,NBR=9 ,NBT=9) 

00040  PARAMETER(NFMX=2,NSGF=2,NB0DMX=19,NSGB=1 ) 

00120  DATA  NN0D0T , NKR0T , NKT0T , IF0RCE/0 ,1,1,1/ 

00130  DATA  GAMA, VISC0,G,T0LR/-0. 78539816, 1.E-6, 9. 80621, 1.E-4/ 

26.  The  PARAMETER  statements  consist  of  the  following  parameters: 

Parameter 

Name  _ Definition _ 

NELE  Total  number  of  elements 

NN0D  Total  number  of  nodes 

NBD  Bandwidth,  which  is  the  maximum  value  of  NBP,  NBT,  and  the  maximum 

value  of  nodal  number  difference  in  each  element.  (The  value  of 
NBD  usually  cannot  be  easily  prepared,  and  an  estimated  value 
should  be  used  for  the  first  run.  The  correct  value  of  NBD  will 
be  in  the  output  of  each  run.  For  example,  in  the  output  file  in 
Appendix  A,  the  correct  value  for  NBD  is  given  Just  before  wave 
period  information.) 

(Continued) 
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Parameter 

Name 


Definition 


NKR 

NKT 

NBR 

NBT 

NFMX 

NSGF 

NBODMX 

NSGB 

27. 

Data  Name 
NNODOT 

NKROT 
NKTOT 
I FORCE 

GAMA 

VISCO 


28. 

01980 
01990  C 
02000  C 
02010  C 
02020 
02030 


Total  number  of  wave  modes  for  reflection  ^k*  in  Equation  14^ 

Total  number  of  wave  modes  for  transmission  ^k~  in  Equation  14^ 

Number  of  nodes  along  the  reflection  boundary  3A+ 

Number  of  nodes  along  the  transmission  boundary  3  A- 
Maximum  nodal  number  allowed  in  a  free  surface  segment  3F 
Number  of  segments  of  free  surface  3F 
Maximum  nodal  number  allowed  in  a  solid  body 
Number  of  bodies  for  force  calculation 

The  DATA  statements  consist  of  the  following  parameters: 

_ Definition _ 

Index  for  the  nodal  value  output 
0  for  no  print  out 

NN0D  for  printing  out  all  the  nodal  potential  values 

Number  of  the  reflection  coefficients  R+  in  Equation  14 
to  be  printed  out 

Number  of  the  transmission  coefficients  R~  in  Equation  14 
to  be  printed  out 

Index  for  force  calculation 

IFORCE=0  for  no  force  calculation  (Otherwise  force  will  be 
calculated. ) 

Y  in  Equation  19  (This  is  the  phase  difference.  If  it  is  not 
known  it  is  usually  taken  to  be  -ir/4  .) 

v  in  Equation  19  (This  is  the  water  viscosity.  If  it  is  not 
known  it  is  usually  taken  to  be  10"^  . ) 

Subroutine  DATAIN 

The  listing  for  subroutine  DATAIN  is  as  follows: 

SUBROUTINE  DATAIN(IG0,0MGSG,0MGA,G) 

0MGSG=(WAVE  FREQUENCY ) **2 /GRAVITY  CONSTANT. 


READ(5,4)  IG0,WAVT 
4  F0RMAT( 1 10, FI 0. 1 ) 
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02040  IF( IGO.EQ.O)  RETURN 

02050  WRITE(6,8)  WAVT 

02060  8  FORMAT( //20( " — " )/ 1 OX , "WAVE  PERIOD,  WAVT=",F8.4," 

02070  1 SECONDS  ,/20('' — ")//) 

02080  OMG A=6 .2831 853/WAVT 

02090  OMGSG=OMGA»OMGA/G 

02100  RETURN 

02110  END 


29.  This  subroutine  requires  the  wave  period  WAVT  be  read  in  with 
unit  of  seconds  as  the  wave  parameter.  If  another  type  of  wave  parameter  is 
chosen  to  be  read  in,  some  modified  statements  to  calculate  OMGA  and  OMGSG  are 
needed  in  the  subroutine.  Also,  the  format  statement  must  be  changed  to  get 
the  right  units.  OMGA  is  the  radian  wave  frequency  u  in  Equation  4  and  is 
2ir/WAVT  ;  OMGSG  =  w2/g  . 

Data  Card  Structure 

30.  Data  are  read  by  the  main  program  and  the  subroutine  DATAIN  from 
device  #5  (TAPE5).  The  following  data  sets  must  all  appear  and  be  placed  in 
order . 


a.  Data  set  1 . 

00180  READ(5,8)  TITLE 

00190  8  F0RMAT(20A4) 

TITLE  has  a  dimension  of  20  and  is  reserved  for  the  title  of 
the  Job.  It  uses  a  single  card  on  which  the  user  types  a  title 
or  identifying  label  anywhere  using  any  or  all  of  columns  1-80. 

b.  Data  set  2. 

00290  READ(5,12)  ( I ,X( I ) ,Z(I ) ,BETA( I ) , J=1 ,NN0D) 

00300  12  F0RMAT(3(I4,3F6.0)) 

In  this  data  set,  node  numbers,  their  x  and  z  coordinates,  and 
friction  coefficients  are  read  in.  I  is  the  nodal  number,  X(I) 
and  Z(I)  are  the  x  and  z  coordinates  of  nodal  I,  respectively, 
and  BETA(I)=  0  is  the  friction  coefficient  as  seen  in  Equa¬ 
tion  19.  Each  card  contains  information  for  three  nodes.  Co¬ 
ordinate  values  must  be  read  in  for  all  nodes.  The  number  of 
cards  read  is  NN0D/3  ,  plus  1  if  there  is  a  remainder. 

c.  Data  set  3. 

00160  2  F0RMAT(20I4) 

00370  READ(5,2)  (J,(IC0N(I,J),I=1,3),L=1,NELE) 
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In  this  data  set,  element  numbers  and  their  nodal  connectivi¬ 
ties  are  read  in.  J  is  the  element  number,  and  ICON(I,J)  are 
the  nodal  connectivities  surrounding  element  J.  Each  card  con¬ 
tains  information  for  five  elements.  Nodal  numbers  must  be 
read  in  counterclockwise  order  for  all  elements.  The  number  of 
cards  read  in  is  NELE/5,  plus  1  if  there  is  a  remainder. 

d.  Data  set  4. 

00160  2  F0RMAT(20I4) 

00610  DO  35  1=1, NSGF 

00620  READ(5,2)  JJ, ( INF( J, I ) ,  J=1,JJ) 

00680  35  CONTINUE 

This  data  set  is  repeated  NSFG  times  to  account  for  NSFG  seg¬ 
ments  on  the  free  surface  3F  .  In  this  data  set,  the  number 
of  nodes  on  each  segment  of  the  free  surface,  along  with  the 
nodal  connectivity  on  the  segment,  are  read  in.  JJ  is  the  num¬ 
ber  of  nodes  on  the  segment  of  the  free  surface,  and  INF(J,I) 
is  the  nodal  connectivity  which  must  be  in  the  order  of  from 
left  to  right. 

e.  Data  set  5. 

00160  2  F0RMAT(20I4) 

00170  4  FORMAT (8F1 0.4) 

00820  DO  45  1=1, NSGB 

00830  READ(5,2)  JJ,(INB0D(J,I),J=1,JJ) 

00840  READ(5,4)  SCL(I),XC(I),ZC(I) 

00910  45  CONTINUE 

If  NSGB=0,  this  data  set  is  skipped.  This  data  set  is  repeated 
NSGB  times  to  account  for  NSGB  segments  on  the  solid  boundaries 
3B  on  which  the  forces  are  to  be  calculated.  On  card  00820, 
the  number  of  nodes  and  the  nodal  connectivities  around  each 
segment  of  the  solid  boundaries  is  read  in.  JJ  is  the  number 
of  nodes  on  the  segment,  and  INB0D(J,I)  is  the  nodal  connectiv¬ 
ity  which  must  be  in  the  sequence  to  make  the  water  domain  area 
positive  (the  water  domain  must  be  on  the  right  hand  side  in 
the  direction  of  nodal  connectivity).  Twenty  values  can  be 
read  from  each  card.  Card  00830  should  be  immediately  followed 
by  card  00840.  On  this  card  the  length  scale  SCL(I),  which  is 
used  to  normalize  the  forces,  and  XC(I)  and  ZC(I),  which  are 
the  reference  center  for  moment  ,  are  read  in.  I  is  the 
segment  of  the  solid  boundary.  y 

f .  Data  set  6. 

01980  SUBROUTINE  DATAIN( IG0,0MGSG,0MGA,G) 

02020  READ(5,4)  IG0,WAVT 

02030  4  F0RMAT( 1 10, F 10. 1 ) 

As  mentioned  in  paragraphs  28  and  29,  this  data  set  consists  of 
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a  set  of  cards  containing  an  integer  number  for  the  sequential 
number  of  wave  periods  and  a  wave  period  in  seconds  indicated 
by  IGO  and  WAVT,  respectively.  IGO  should  never  be  equal  to 
0  . 

g.  Data  set  7. 

At  least  one  card  with  0  typed  in  column  10  should  be  read  in 
to  indicate  the  end  of  all  the  input  data. 

31.  All  printed  output  is  directed  to  device  #6  (tape  6). 
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PART  III:  EXAMPLE  PROBLEM 


32.  As  an  example  a  fixed,  semi- immersed  circular  cylinder  in  a  finite 
water  depth  (as  shown  in  Figure  4),  along  with  the  finite  element  mesh,  will 
be  used  here.  For  clarity,  the  length  unit  in  Figures  4  through  6  is  replaced 
by  metres  instead  of  inches  in  this  calculation.  In  the  water  domain  (A), 
there  are  81  nodes  and  116  elements.  There  exists  one  fixed  floating  ob¬ 
stacle.  There  are  two  free  surfaces,  one  along  nodes  9  and  20  and  the  other 
along  nodes  72  and  81.  To  run  this  the  program,  the  values  in  the  two 
PARAMETER  and  two  DATA  statements  in  the  main  program  (given  in  paragraph  33) 
need  be  changed  along  with  the  input  data  in  paragraph  34.  The  output  data 
are  given  in  Appendix  A. 


PARAMETER  and  DATA  Statements 


33.  In  the  main  program  the  necessary  changes  needed  to  implement  this 
particular  problem  are  found  in  the  PARAMETER  statements  and  the  DATA  state¬ 
ments.  The  values  are  as  follows: 

00030  PARAMETER(NELE=116,NN0D=81,NBD=13,NKR=6,NKT=6,NBR=9,NBT=9) 

00040  PARAMETER( NFMX=2 , NSGF=2 , NB0DMX= 1 9 , NSGB=1 ) 

00120  DATA  NN0D0T,NKR0T,NKT0T, I F0RCE/8 1,6, 6,1/ 

00130  DATA  GAMA, VISC0,G,T0LR/-0. 7853981 6, 1.E-6, 9. 8062 1,1.E-4/ 

Input  Data  File 

34.  The  input  file  is  given  in  the  following  list  with  accompanying 
card  numbers.  The  title  is  on  card  00010.  All  subsequent  cards  are  assumed 
to  have  exactly  80  columns.  Cards  00020  through  00280  represent  the  nodes  and 
their  x  and  z  coordinates  and  friction  coefficients.  The  elements  and  their 
nodal  connectivities  are  on  cards  00290  through  00520.  Cards  00530  and  00540 
have  the  number  of  nodes  and  the  nodal  connectivities  on  each  segment  of  free 
surface.  Card  00550  has  the  number  of  nodes  and  the  nodal  connectivities  on 
each  segment  of  the  solid  body  for  force  calculation.  Card  00560  has  the 
length  scale  and  the  reference  center  chosen  to  be  the  radius  and  the  center 
of  the  semicircular  cylinder.  Cards  00570  and  00580  have  the  index  and  the 
wave  period.  Card  00590  has  a  zero  to  represent  the  end  of  input. 
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Input  File  Listing 


00010 

REFLECTION.  TRANSMISSION* 

AND 

NAVE  FORCES  FOR  LINEAR  WAVE 

00020 

L 

-8.0  -28.0 

0.00  2  - 

•8.0 

-22.0  0.00 

3 

-8.0  -17.0 

0.00 

00030 

4 

-8.0  -13.0 

0.00  5  ■ 

-8.0 

-10.0  0.00 

6 

-8.0  -7.5 

0.00 

00040 

7 

-8.0  -5.0 

0.00  8  ■ 

•8.0 

-2.5  0.00 

9 

-8.0  0.0 

0.00 

00050 

10 

-6.0  -28.0 

10.00  11  ■ 

•6.0 

-22.0  0.00 

12 

-6.0  -17.0 

0.00 

00060 

13 

-6.0  -13.0 

0.00  14  ■ 

-6.0 

-10.0  0.00 

15 

-6.0  -7.5 

0.00 

00070 

16 

-6.0  -5.0 

0.00  17-5, 

,196- 

■3.000  10.00 

18-5.638-2.052 

10,00 

00080 

19- 

■5.909-1.042 

10.00  20  - 

■6.0 

0.0  10.00 

21 

-4.0  -28.0 

10. 00 

00090 

')'> 

-4,0  -22.0 

0.00  23  ■ 

-4.0 

-17.0  0.00 

24 

-4.0  -13.0 

o.oc 

00100 

25 

-4.0  -10.0 

0.00  26  - 

■4.0 

-7.5  0.00 

27-3.857-4.596 

10.00 

00110 

28- 

-4.596-3.857 

10.00  29  - 

■2.0 

-28.0  10.00 

30 

-2.0  -22.0 

0.00 

00120 

31 

-2.0  -17.0 

0.00  32  • 

-2.0 

-13.0  0.00 

33 

-2.0  -10.0 

0.00 

00130 

34 

-2.0  -7.5 

0.00  35-3, 

,000- 

■5.196  10.00 

36- 

■2.052-5.638 

10.00 

00140 

37- 

•1.042-5.909 

10.00  38 

0.0 

-28.0  10.00 

39 

0.0  -22.0 

0,00 

00150 

40 

0.0  -17.0 

0.00  41 

0.0 

-13.0  0.00 

42 

0.0  -10.0 

0.00 

00160 

43 

0.0  -7.5 

0.00  44 

0.0 

-6.0  10.00 

45 

2.0  -28.0 

10.00 

00170 

46 

2.0  -22.0 

0.00  47 

2.0 

-17.0  0.00 

48 

2.0  -13.0 

0.00 

00180 

49 

2.0  -10.0 

0.00  50 

2.0 

-7.5  0.00 

51 

1.042-5.909 

10. CO 

00190 

52 

2.052-5.638 

10.00  53  3.000- 

■5.196  10.00 

54 

4.0  -28.0 

10.00 

00200 

55 

4.0  -22.0 

0.00  56 

4.0 

-17.0  0.00 

57 

4.0  -13.0 

0.00 

00210 

58 

4.0  -10.0 

0.00  59 

4.0 

-7.5  0,00 

60 

3.857-4.596 

10.00 

00220 

61 

4.596-3.857 

10.00  62 

6.0 

-28.0  10.00 

63 

6.0  -22.0 

0.00 

00230 

64 

6.0  -17.0 

0.00  65 

6.0 

-13.0  0.00 

66 

6.0  -10.0 

0.00 

00240 

67 

6.0  -7.5 

0.00  68 

6.0 

-5.0  0.00 

69 

5.196-3.000 

10.00 

00250 

70 

5.638-2.052 

10.00  71  5. 

909- 

1.042  10.00 

72 

6.0  0.0 

10.00 

00260 

73 

8.0  -28.0 

0.00  74 

8.0 

-22.0  0.00 

75 

8.0  -17.0 

0.00 

00270 

76 

8.0  -13.0 

0.00  77 

8.0 

-10.0  0.00 

78 

8.0  -7.5 

0.00 

00280 

79 

8.0  -5.0 

0.00  80 

8.0 

-2.5  0.00 

81 

8.0  0.0 

0.00 

00290 

1 

1  10  2 

2  2  10 

11 

3  2  11 

3 

4  3  11 

12  5 

3 

12 

4 

00300 

6 

4  12  13 

7  4  13 

5 

8  5  13 

14 

9  5  14 

6  10 

6 

14 

15 

00310 

11 

6  15  7 

12  7  15 

16 

13  7  16 

8 

14  8  16 

17  15 

3 

17 

n 

00320 

16 

3  18  19 

17  8  19 

9 

18  9  1® 

20 

19  10  21 

11  20 

11 

21 

00330 

21 

11  22  12 

22  12  22 

23 

23  12  23 

13 

24  13  23 

24  25 

13 

24 

14 

00340 

26 

14  24  25 

27  14  25 

15 

28  15  25 

26 

29  15  26 

■6  30 

16 

.26 

A.  ' 

00350 

31 

16  27  28 

32  16  28 

17 

33  21  29 

22 

34  22  29 

30  35 

22 

30 

23 

00360 

36 

23  30  31 

37  23  31 

24 

38  24  31 

32 

39  24  32 

25  40 

25 

32 

33 

00370 

41 

25  33  26 

42  26  33 

34 

43  26  34 

35 

44  26  25 

27  45 

34 

56 

35 

00380 

46 

29  38  30 

47  30  38 

39 

48  30  39 

31 

49  31  39 

40  50 

31 

40 

32 

00390 

51 

32  40  41 

52  32  41 

33 

53  33  41 

42 

54  33  42 

34  55 

34 

42 

43 

00400 

56 

34  43  37 

57  34  37 

36 

58  43  44 

37 

59  38  45 

46  60 

38 

46 

39 

00410 

61 

39  46  47 

62  40  39 

47 

63  40  47 

48 

64  41  40 

48  65 

41 

48 

4° 

00420 

66 

42  41  49 

67  42  49 

50 

68  43  42 

50 

69  43  50 

51  70 

43 

51 

44 

00430 

71 

51  50  52 

72  45  54 

55 

73  45  35 

46 

74  46  55 

56  75 

47 

46 

56 

00440 

76 

47  56  57 

77  48  47 

57 

78  48  57 

58 

79  49  48 

58  80 

49 

58 

59 

00450 

81 

50  49  59 

82  50  59 

53 

83  52  50 

53 

84  53  59 

60  85 

54 

62 

63 

00460 

86 

55  54  63 

87  55  63 

64 

88  56  55 

64 

89  56  64 

65  90 

57 

56 

65 

00470 

91 

57  65  66 

92  58  57 

66 

93  58  66 

67 

94  59  58 

67  95 

59 

67 

63 

00480 

96 

60  59  68 

97  60  68 

61 

98  61  68 

69 

99  62  73 

74  100 

63 

62 

74 

00490 

101 

63  74  75 

102  64  63 

75 

103  64  75 

76 

104  65  64 

76  105 

65 

76 

“*7 

00500 

106 

66  65  77 

107  66  77 

78 

108  67  66 

78 

109  67  78 

79  110 

63 

67 

79 

00510 

111 

68  79  80 

112  68  80 

69 

113  69  80 

70 

114  70  80 

71  115 

71 

80 

31 

00520 

116 

72  71  81 

00530 

2 

9  20 

00540 

2 

72  81 

00550 

19 

20  19  18 

17  28  27 

35 

36  37  44 

51 

52  53  60 

61  69 

70 

71 

72 

00560  6.0  0.0  0.0 

00570  1  12.0 

00580  2  6.0 

00590  0 
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APPENDIX  A:  OUTPUT  DATA  FILE 


The  program  with  the  correct  PARAMETER  and  DATA  statements  and  the  input 
data  file  was  submitted  and  run  on  a  Control  Data  Cyber  170  Model  760  compu¬ 
ter.  Results  in  the  output  file  are  shown  in  the  following  pages  of  this  ap¬ 
pendix.  In  the  output,  nodal  potential  is  given  in  the  form  of  */(iga0/«)  , 
where  igao/w  is  the  factor  of  the  incident  wave  as  given  in  Equation  13. 

The  reflection  and  transmission  coefficients  are  Rm  and  R*  in  Equation  14, 
and  the  wave  forces  and  moment  coefficients  are  fx  ,  fz  ,  and  niy  as  given 
in  Equations  21,  22,  and  23,  respectively.  The  numerical  results  are  self- 
explanatory.  The  CPU  time  for  this  run  is  less  than  2  sec.  The  output  file 
listing  is  as  follows: 


A1 


REFLECTION*  TRANSMISSION,  AND  WAVE  FORCES  FOR  LINEAR  WAVE 


TOTAL  NUMBER  OF  ELEMENTS,  NELE-  116 
TOTAL  NUMBER  OF  NODES,  NNOD=  81 
TOTAL  NUMBER  OF  EI6EN  VALUES  FOR  REFLECTION,  NKR=  6 
TOTAL  NUMBER  OF  EIGEN  VALUES  FOR  TRANSMISSION,  NKT=  6 
MAX  NODAL  NO.  ALLOWED  IN  A  SEGMENT,  NFMX=  2 
INDEX  FOR  FORCE  CALCULATION  <N0N£=0>»  IFORCE=  1 
MAX  NODAL  NO.  ALLOWED  IN  A  BODY,  NBODMX-  1? 


(X, Z)  COORDINATE  AND  FRICTION  COFF  FOR  EACH  NOD 


NOD 

X 

Z 

FRIC 

NOD 

X 

Z 

FRIC 

NOD 

FRIC 

1 

-8.0 

-28.0 

.0 

2 

-8.0 

-22.0 

.0 

3 

-8.0 

i 

■'J 

o 

.0 

4 

-8.0 

-13.0 

.0 

5 

-8.0 

-10.0 

.0 

<s 

-8.0 

-7.5 

.  0 

7 

-8.0 

-5.0 

.0 

8 

-8.0 

-2.5 

.0 

9 

-3.0 

.0 

.0 

10 

-6.0 

-28.0 

10,0 

11 

-6.0 

-22.0 

.0 

12 

-6.0 

-17,0 

,0 

13 

-6.0 

-13.0 

.0 

14 

-6.0 

-10.0 

.0 

15 

-6.0 

-7.5 

,0 

16 

-6.0 

-5.0 

.0 

17 

-5,2 

-3.0 

10.0 

18 

-5.6 

-2.1 

10.0 

19 

-5.9 

-1.0 

10.0 

20 

-6.0 

,0 

10.0 

21 

-4.0 

-2S.0 

10,0 

22 

-4.0 

-22.0 

.0 

23 

-4.0 

-17,0 

.0 

24 

-4.0 

-13.0 

.0 

25 

-4.0 

-10.0 

.0 

26 

-4,0 

-7,5 

.0 

27 

-3.9 

-4.6 

10,0 

28 

-4.6 

-3.9 

10.0 

29 

-2.0 

-28.0 

10.0 

30 

-2,0 

-22.0 

.0 

31 

-2.0 

-17.0 

.0 

32 

-2.0 

-13.0 

.0 

33 

-2.0 

-10.0 

.0 

34 

-2.0 

-7 . 5 

A 

♦  \J 

35 

-3.0 

-5.2 

10.0 

36 

-2.1 

-5.6 

10.0 

37 

-1.0 

-5.9 

10.0 

38 

.0 

-28.0 

10,0 

39 

.0 

-22.0 

.0 

40 

.0 

-17.0 

.0 

41 

.0 

-13.0 

,0 

42 

,0 

-10.0 

,0 

43 

.0 

-7,5 

.0 

44 

.0 

-6.0 

10,0 

45 

2.0 

-28.0 

10.0 

46 

2.0 

-22.0 

.0 

47 

2.0 

-17.0 

.0 

48 

2.0 

-13.0 

.0 

4V 

2.0 

-10.0 

,0 

50 

2,0 

-7,5 

.0 

51 

1.0 

-5.9 

10.0 

52 

2.1 

-5.6 

10.0 

53 

3.0 

-5.2 

10.0 

54 

4.0 

-28.0 

10,0 

55 

4.0 

-22.0 

.0 

56 

4.0 

-17,0 

.0 

57 

4.0 

-13.0 

.0 

58 

4.0 

-10.0 

.0 

59 

4,0 

-7.5 

,0 

60 

3.9 

-4.6 

10,0 

61 

4.6 

-3.9 

10.0 

62 

6.0 

-28.0 

10.0 

63 

6.0 

-22.0 

.0 

64 

6.0 

-17.0 

,0 

65 

6.0 

-13,0 

.0 

66 

6.0 

-10,0 

.0 

67 

6.0 

-7.5 

.0 

68 

6,0 

-5.0 

,0 

69 

5,2 

-3.0 

10.0 

70 

5.6 

-2.1 

10.0 

71 

5.9 

-1.0 

10.0 

72 

6.0 

.0 

10,0 

73 

8.0 

-28.0 

.0 

74 

8.0 

-22.0 

.0 

75 

8.0 

-17.0 

.0 

76 

8.0 

-13,0 

,0 

77 

8.0 

-10.0 

.0 

78 

8.0 

-7,5 

,0 

79 

8.0 

-5.0 

,0 

80 

8.0 

-2,5 

.0 

81 

8.0 

.0 

.0 

A2 


o  o  o  o 


NODAL  CONNECTIVITY 


ELEM 

N1 

N2 

N3 

ELEM 

N1 

N2 

N3  ELEM 

N1 

N2 

N3 

ELEM 

N1 

N2 

N3 

1 

1 

10 

2 

2 

2 

10 

11 

3 

2 

11 

3 

4 

3 

11 

12 

5 

3 

12 

4 

6 

4 

12 

13 

7 

4 

13 

5 

8 

5 

13 

14 

9 

5 

14 

6 

10 

6 

14 

15 

11 

6 

15 

7 

12 

7 

15 

16 

13 

7 

16 

8 

14 

8 

16 

17 

15 

8 

17 

18 

16 

8 

18 

19 

17 

8 

19 

9 

18 

9 

19 

20 

19 

10 

21 

11 

20 

11 

21 

22 

21 

11 

22 

12 

22 

12 

22 

23 

23 

12 

23 

13 

24 

13 

23 

24 

25 

13 

24 

14 

26 

14 

24 

25 

27 

14 

25 

15 

28 

15 

25 

26 

29 

15 

26 

16 

30 

16 

26 

27 

31 

16 

27 

28 

32 

16 

28 

17 

33 

21 

29 

22 

34 

22 

29 

30 

35 

22 

30 

23 

36 

23 

30 

31 

37 

23 

31 

24 

38 

24 

31 

32 

39 

24 

32 

25 

40 

25 

32 

33 

41 

25 

33 

26 

42 

26 

33 

34 

43 

26 

34 

35 

44 

26 

35 

27 

45 

34 

36 

35 

46 

29 

38 

30 

47 

30 

38 

39 

48 

30 

39 

31 

49 

31 

39 

40 

50 

31 

40 

32 

51 

32 

40 

41 

52 

32 

41 

33 

53 

33 

41 

42 

54 

33 

42 

34 

55 

34 

42 

43 

56 

34 

43 

37 

57 

34 

37 

36 

58 

43 

44 

37 

59 

38 

45 

46 

60 

38 

46 

39 

61 

39 

46 

47 

62 

40 

39 

47 

63 

40 

47 

48 

64 

41 

40 

48 

65 

41 

48 

49 

66 

42 

41 

49 

67 

42 

49 

50 

68 

43 

42 

50 

69 

43 

50 

51 

70 

43 

51 

44 

71 

51 

50 

52 

72 

45 

54 

55 

73 

45 

55 

46 

74 

46 

55 

56 

75 

47 

46 

56 

76 

47 

56 

57 

77 

48 

47 

57 

78 

48 

57 

58 

79 

49 

48 

58 

80 

49 

58 

59 

81 

50 

49 

59 

82 

50 

59 

53 

83 

52 

50 

53 

84 

53 

59 

60 

85 

54 

62 

63 

86 

55 

54 

63 

87 

55 

63 

64 

88 

56 

55 

64 

89 

56 

64 

65 

90 

57 

56 

65 

91 

57 

65 

66 

92 

58 

57 

66 

93 

58 

66 

67 

94 

59 

58 

67 

95 

59 

67 

68 

96 

60 

59 

68 

97 

60 

68 

61 

98 

61 

68 

69 

99 

62 

73 

74 

100 

63 

62 

74 

101 

63 

74 

75 

102 

64 

63 

75 

103 

64 

75 

76 

104 

65 

64 

76 

105 

65 

76 

77 

106 

66 

65 

77 

107 

66 

77 

78 

108 

67 

66 

78 

109 

67 

78 

79 

110 

68 

67 

79 

111 

68 

79 

80 

112 

68 

80 

6'.’ 

113 

69 

80 

70 

114 

70 

80 

71 

115 

71 

80 

81 

116 

72 

71 

81 

NUMBER  OF  NODES  ON  REFLECTION  DOMAIN?  NBR=  9 
THEIR  CONNECTIVITY  ARE! 

73  74  75  76  77  78  79  80  81 


NUMBER  OF  NODES  ON  TRANSMISSION  DOMAIN?  NBT=  9 
THEIR  CONNECTIVITY  AREt 
123456789 


9  20 


NUMBER  OF  SEGMENTS  OF  FREE  SURFACE?  NSGF-  2 

NUMBER  OF  NODES  ON  1-TH  SEGMENT?  NF(  1)-  2 

THEIR  CONNECTIVITY  ARE  I 


NUMBER  OF  NODES  ON  2-TH  SEGMENT?  NF(  2)-  2 

THEIR  CONNECTIVITY  ARE? 


A3 


72  81 


FOR  REFLECTION  DOMAIN i 

WATER  DEPTH.  HR=  28,00  HORIZONTAL  EXTENT,  XR=  8.00 
FOR  TRANSMISSION  DOMAIN,* 

WATER  DEPTH,  HT=  28,00  HORIZONTAL  EXTENT,  XT=  -8,00 


NUMBER  OF  BODY  FOR  FORCE  CALCULATION,  NSGB^ 


NUMBER  OF  NODES  ON  1-TH  BODY,  NBOD <  1)  =  19 

LENGTH  SCALE,  SCL(I)=  6.00 

REFERENCE  CENTER,  <XC,ZC)=  .00 

THEIR  CONNECTIVITY  ARE! 

20  19  18  17  28  27  35  36  37  44  51  52  53  60  61 


,00 


69  70  71 


BANDWITH,  NBD=  13 


WAVE  PERIOD,  WAVT-  12.00  SECONDS 


-  THE  SOLUTION  OF  THE  SYSTEM  -  - 


NODAL-POTENTIAL 


NODE 

REAL-PART 

IMAGE-PART 

ABS-VALUE 

PHASE 

1 

.6329 

.0459 

,6346 

,072 

2 

.6428 

.0504 

.6448 

.078 

3 

,6660 

.0644 

.6691 

.096 

4 

.6923 

.0842 

.6984 

.121 

5 

.7197 

,1098 

.7280 

,151 

6 

.7461 

.1390 

.7589 

.184 

AH 


w 


7 

.7775 

,1767 

.7974 

.223 

8 

.8202 

.2212 

.8495 

.263 

9 

.8762 

.2532 

,9121 

.281 

10 

.6459 

-.0053 

.6459 

-.008 

11 

.6557 

-.0022 

.6557 

-.003 

12 

.6780 

.0080 

.6781 

.012 

13 

.7040 

.0243 

.7044 

.034 

14 

.7276 

.0467 

.7291 

.064 

15 

.7499 

.0766 

.7538 

.102 

16 

.7756 

.1235 

.7853 

.  158 

17 

.8011 

.1544 

.8158 

.190 

18 

.8215 

.1890 

.8429 

.226 

19 

,8462 

.2143 

.8729 

.248 

20 

.8719 

.2286 

,9014 

.256 

21 

.6569 

-.0575 

.6594 

-.087 

22 

,6665 

-.0559 

.6689 

-.084 

23 

.6883 

-.0503 

.6902 

-.073 

24 

.7130 

-.0396 

.7141 

-.055 

25 

.7341 

-.0228 

.7344 

-.031 

26 

.7514 

.0032 

.7515 

.004 

27 

.7712 

.0593 

.7735 

.077 

28 

.7832 

.1110 

.7910 

.141 

29 

.6657 

-.1099 

.6747 

-.164 

30 

.6753 

-.1102 

.6843 

-.162 

31 

.6970 

-.1098 

.7056 

-.156 

32 

,7212 

-.1068 

.7290 

-.147 

33 

.7406 

-.0993 

.7472 

-.133 

34 

,7545 

-.0844 

.7592 

-.111 

35 

,7637 

.0055 

.7638 

.007 

36 

.7607 

-.0532 

.7626 

-.070 

37 

.7616 

-.1188 

.7708 

-.155 

38 

.6724 

-.1624 

.6918 

-.237 

39 

.6822 

-.1647 

.7017 

-.237 

40 

.7042 

-.1700 

,7244 

-.237 

41 

.7288 

-.1760 

.7497 

-.237 

42 

.7485 

-.1808 

.7701 

-.237 

43 

.7624 

-.1841 

.7843 

-.237 

44 

.7665 

-.1851 

.7886 

-.237 

45 

.6770 

-.2142 

.7101 

-.306 

46 

.6869 

-.2186 

.7209 

-.308 

47 

.7098 

-.2298 

.7460 

-.313 

48 

.7358 

-.2450 

.7756 

-.321 

49 

.7580 

-.2626 

,8022 

-.334 

50 

.7758 

-.2852 

,8266 

-.352 

51 

.7757 

-.2525 

,8158 

-.315 

52 

.7892 

-.3211 

.8520 

-.386 

53 

.8052 

-.3844 

.8923 

-.445 

54 

.6793 

-.2652 

.7292 

-.372 

55 

.6896 

-.2715 

.7411 

-.375 

56 

.7136 

-.2882 

.7696 

-.384 

57 

.7419 

-.3117 

.8048 

-.398 

58 

.7678 

-.3399 

.8396 

-.417 

59 

.7918 

-.3759 

.8765 

-.443 

60 

.8250 

-.4447 

.9372 

-.494 

A5 


61 

.8489 

-.5050 

.9877 

-.537 

62 

.6794 

-.3147 

.7487 

-.434 

63 

.6900 

-.3227 

.7617 

-.437 

64 

.7155 

-.3444 

.7941 

-.449 

65 

.7463 

-.3745 

.8350 

-.465 

66 

.7761 

-.4098 

.8777 

-.486 

67 

.8061 

-.4522 

.9243 

-.511 

68 

.8435 

-.5144 

.9880 

-.548 

69 

.8772 

-.5594 

1.0404 

-.568 

70 

.9064 

-.6058 

1.0902 

-.589 

71 

.9380 

-.6445 

1.1381 

-.602 

72 

,9684 

-.6721 

1.1787 

-.607 

73 

.6769 

-.3623 

.7677 

-.491 

74 

.6879 

-.3717 

.7819 

-.495 

75 

.7150 

-.3978 

.8182 

-.508 

76 

.7481 

-.4323 

.8640 

-.524 

77 

.7815 

-.4722 

.9131 

-.544 

78 

.8157 

-.5160 

.9652 

-,5M 

79 

.8571 

-.5713 

1 .0300 

-.588 

80 

.9119 

-.6392 

1.1136 

-.611 

81 

.9781 

-.7003 

1.2030 

-.621 

COEFFICIENTS 

FOR  REFLECTION  DOMAIN 

COEF 

REAL-PART 

IMAGE-PART 

ABS-VALUE  PHASE 

1 

-.0436 

-.3319 

.3347  -1.701 

2 

-.0766 

-.0995 

,1255  -2.227 

3 

-.0144 

-.1196 

,1204  -1.691 

4 

,0088 

-.1285 

.1288  -1.50 > 

5 

.0302 

-.1431 

•  1463  - 1  ♦ 

6 

.0530 

-.1554 

♦1642  -1*242 

COEF 

COEFFICIENTS 

REAL-PART 

FOR  TRANSMISSION 
IMAGE-PART  ABS- 

DOMAIN 

-VALUE 

PHASE 

1 

.9334 

-.1243 

.9416 

-.132 

2 

-.1021 

.1429 

.1757 

2.191 

3 

-.0430 

.1338 

.1405 

1.882 

4 

-.0191 

.1321 

.1335 

1.714 

5 

-.0018 

.1376 

.1377 

1.584 

6 

,0197 

.1415 

.1429 

1.432 

WAVE  FORCES 

J  INDX»1=FX» 

2=FZ»3=MY 

INDX 

REAL-PART 

IMAGE-PART 

A8S-VALU E 

PHASE 

1 

-.0711 

.6551 

.6589 

1.679 

2 

1.5950 

-.3855 

1.6409 

-.237 

3 

-.0014 

-.0009 

.0017 

-2.566 
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WAVE  PERIOD »  WAVT= 


6.00  SECONDS 


-  THE  SOLUTION  OF  THE  SYSTEM  -  - 


NODAL- POTENTIAL 


NODE 

REAL-PART 

IMAGE-PART 

ABS-VALUE 

PHASE 

1 

.1117 

-.0860 

.1410 

-.656 

.1171 

-.0828 

.1434 

-.615 

3 

.1288 

-.0682 

.1457 

-.487 

4 

.1420 

-.0418 

.1481 

-.286 

5 

.1536 

-.0023 

.1536 

-.015 

6 

.1653 

.0474 

.1720 

.279 

7 

.1811 

.1161 

.2151 

.570 

8 

.2105 

.2070 

.2952 

.7/7 

9 

.2675 

.2999 

.4018 

.842 

10 

.1244 

-.1043 

.1623 

-.698 

11 

.1314 

-.1042 

.1677 

-.670 

12 

.1469 

-.0985 

.1769 

-.591 

13 

.1643 

-.0824 

.1838 

-.465 

14 

.1780 

-.0532 

.1858 

-.290 

15 

.1887 

-.0098 

.1889 

-,o:,2 

16 

.1983 

.0610 

,2074 

.298 

17 

.2157 

.1187 

.2462 

.503 

18 

.2268 

.1754 

.2868 

.658 

19 

.2477 

.2271 

.3361 

.742 

20 

.2757 

.2681 

.3646 

.772 

21 

.1356 

-.1228 

.1829 

-.736 

22 

,1442 

-.1264 

.1918 

-.720 

23 

.1640 

-.1314 

,2102 

-.675 

24 

.1862 

-.1296 

.2268 

-.608 

25 

.2029 

-.1151 

.2333 

-.516 

26 

.2128 

-.0846 

,2290 

-.378 

27 

.2149 

-.0110 

.2152 

-.051 

28 

.2108 

.0553 

.2179 

,257 

29 

.1448 

-.1404 

,2017 

-.7/0 

30 

.1552 

-.1483 

,2147 

-.763 

31 

.1797 

-.1655 

.2443 

-.744 

32 

.2080 

-.1817 

.2762 

-.718 

33 

.2302 

-.1878 

.2971 

-.684 

34 

.2432 

-.1792 

.3021 

-.635 

A7 


35 

.2246 

-.0758 

.2370 

-.325 

36 

.2404 

-.1438 

.2801 

-.539 

37 

.2630 

-.2204 

.3432 

-.697 

33 

.1517 

-.1562 

.2178 

-.800 

39 

.1639 

-.1687 

.2353 

-.800 

40 

.1935 

-.1992 

. 2777 

-.800 

41 

.2294 

-.2362 

,3293 

-.800 

42 

.2608 

-.2685 

.3743 

-.800 

43 

.2842 

-.2926 

.4079 

-.800 

44 

.2914 

-.3000 

.4182 

-.800 

45 

.1559 

-.1693 

.2301 

-.826 

46 

.1698 

-.1863 

.2521 

-.831 

47 

.2045 

-.2299 

.3077 

-.844 

48 

.2492 

-.2889 

.3816 

-.859 

49 

.2929 

-.3507 

.4570 

-.875 

50 

.3339 

-.4150 

.5326 

-.893 

51 

.3272 

-.3873 

,5070 

-.869 

52 

.3725 

-.4870 

,6131 

-.918 

53 

.4227 

-.5903 

.7261 

-.949 

54 

.1570 

-.1785 

.2378 

-.849 

55 

.1723 

-.1995 

.2636 

-.858 

56 

.2117 

-.2554 

.3317 

-.8:-’9 

57 

.2651 

-.3350 

.4272 

-.901 

58 

.3224 

-.4256 

.5339 

-.922 

59 

,3841 

-.5298 

.6544 

-.943 

60 

.4832 

-.7073 

.8566 

-.971 

61 

.5581 

-.8461 

1.0136 

-.988 

62 

.1547 

-.1832 

.2398 

-.869 

63 

.1709 

-.2071 

.2685 

- .  8;>  1 

64 

.2141 

-.2731 

.3470 

-.906 

65 

.2745 

-.3692 

.4601 

-.931 

66 

.3437 

-.4837 

.5934 

-.953 

67 

.4238 

-.6203 

.7513 

-.971 

68 

.5363 

-.8167 

,9770 

-.990 

69 

.6506 

-1,0093 

1.2009 

-.998 

70 

.7491 

-1.1782 

1.3962 

-1,001 

71 

.8642 

-1.3686 

1,6186 

-1.008 

72 

.9821 

-1.5582 

1.8419 

-1.008 

73 

.1489 

-.1825 

.2355 

-.887 

74 

.1652 

-.2080 

.2656 

-.900 

75 

.2107 

-.2812 

.3513 

-.928 

76 

.2748 

-.3873 

,4749 

-.954 

77 

.3519 

-.5178 

.6261 

-.974 

78 

.4426 

-.6727 

.8052 

-.989 

79 

.5668 

-.8852 

1,0511 

-1.001 

80 

.7511 

-1.1955 

1.4119 

-1,010 

81 

1.0029 

-1.6039 

1.8916 

-1.012 

COEFFICIENTS 

FOR  REFLECTION  DOMAIN 

COEF 

REAL-PART 

IMAGE-PART 

AES-VALUE 

PHASE 

1 

-.3837 

-.8337 

.9177 

-2.002 

A8 


2 

-.0699 

.0454 

.0834  2.565 

3 

-.0084 

-.0541 

.0548  -1.725 

4 

.0416 

-.1225 

.1294  -1,245 

5 

.0898 

-.1914 

.2114  -1,132 

6 

.1470 

-.2702 

.3076  -1.073 

COEF 

COEFFICIENTS 

REAL-PART 

FOR  TRANSMISSION  DOMAIN 
IMAGE-PART  ABS-VALUE 

PHASE 

1 

.3561 

-.1649 

.3925 

-.434 

2 

-.1047 

.1348 

.1707 

2.231 

3 

-.0918 

.1585 

.1832 

2.096 

4 

-.0660 

.1508 

.1646 

1.984 

5 

-.0474 

.1528 

.1600 

1.872 

6 

-.0278 

.1590 

O- 

*-* 

at 

1.744 

WAVE  forces:  INDX,1=FX 

>2=FZ»3=HY 

INDX 

REAL-PART 

IMAGE-PART 

ABS-VALUE 

PHASE 

i 

-.4223 

1.0905 

1.1694 

1.940 

2 

.7025 

-.7240 

1.0088 

-.800 

3 

-.0012 

-.0033 

.0035 

-1.926 

****»****cpu  time,  dt=  1.9402E+00  seconds  ********** 
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APPENDIX  B:  WAVERTF  PROGRAM  LISTING 


The  program  WAVERTF  (WAVE  Reflection,  Transmission,  and  Forces)  is  writ¬ 
ten  in  FORTRAN  Extended  IV  or  FORTRAN  V.  All  FORTRAN  callable  subroutines  are 
self-contained  in  the  program.  FORTRAN  callable  functions  that  are  used  are 
only  those  that  the  system  recognizes.  The  program  is  listed  on  the  following 
pages . 
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00010  PROGRAM  WAVER  F  •;  INPUT  “OUTPUT  ?  TAPES-  IHFVT  ?  TAP Ec -OUTPUT) 

00020  C  WES-CERC  PROGRAM.  HSCHEN.  1984 

00030  PARAMETER ( NELE~1 16  ?  NN0P=81  ?  NBD=  13?  NKR=6 * NKT =6  ?  NBRa? ? NBT-9  ’■ 

00040  PARAMETER < NFMX-2 • N3GF  =2 ? NB0DMX*1? - MS38= 1 • 

00 )50  COMPLEX  3YSK < NNOD ?  NBD > • SYSO ( NNOD 5 ? SYSBR < NKF ) • 3 YSKR ' NKR •  ti£R ) : 

00060  1SYS0R ( NKR >  ?  SYSD T < NKT  i ? SYSKT ( NKT,NBR>? SYS3T { NKT > .» EPSN < NNOD .  ? 

00070  2XIAMD<NNQD) 

00080  DIMENSION  11TLE!20> » ICON ( 3 >  NELE )  ?  X ! NNOD >  * Z < NNOD )  ?  INBRO'BR* » 

00990  !  I  NUT  { NUT ) ? NF  <  N8GF ) » INF  •  NFMX ? NS6F  > ,  WKR  ( NKR ) ? W  KR  ( NKR ) ?  WKT  ■■  NX’ 

00100  2WKHT(NKT)fSCL < NSGB ) ? XC ( NSGB > ?ZC(NSGD) ?MBODl NS5B) ? 

00 1 1 0  3  INBOD  •:  NBODMX  >  NSGB  '■  >  BE  T  A !  NNOD ) 

001 20  DATA  NNODOT  ?  NKROT ? NKTOT >  IF0RCE/B1 ? 3 .•  6>  1  / 

00170  DATA  G AM A  ?  V I SCO  ?  G  ?  T9LR /-0 « 78539B 1 6  ? 1 • E- 6  ?  9 .  8062 1 ? l.E-4/ 

00140  CC  CALL  CRT IMF • CA) 

00150  PI-3. 1415926 

00160  2  FORMAT < 2014 ) 

OOl'A-  4  F0RMAf<8Fll>.4 ') 


00180 

RF  AD  ( 5  ?  8 )  i  I  'M.F 

00190 

3  FORMA T(20A4> 

00200 

WRITE! 6: .9)  Dll  F 

00710 

9  FORMAT ■ //26t ' . 

00220 

WKITE(6?101  NELE 

00230 

0  0  2  *i  0 
r«  •"}  c  /i 

00260 

r.rviic, 


J/2CA4/2 6C-  -*)) 

NNOD • NKR ? NK T * NFMX » I FORCE ? NBODMX 
10  FORMAT <//24X*’TGTAL  NUMBER  OF  ELEMENTS?  NELE**? 14/ 

1  27X? ' TOTAL  NUMBER  OF  NODES?  NNOD- '• H/'V. "TOTAL  NUMBER  Or  EIGEN  0 
2ALUES  FOR  REFLECTION?  NKR=" ? I4/4X. *  TOTAL  NUMBER  OF  EI3EM  MA.UEB  :  C 
3K  TRANSMISSION*  NKT - ’ ? 1 4/1 4x  ?  *MA'<  NODAL  NO.  ALLOWED  IN  A  3ESMEN • * 
4NFMX“*»I4/iOy ?* INDEX  FOR  FORCE  CALCULATION  (N9NE=0).  -FORCE-*  ■•I 


00230 

5  15X? ’MAX  NODAL  NO.  ALLOWED  IN  A  BODY*  NBODMX- * ? 

1 4  ;  • 

002  90 

REAP <5?  12.-  <  1  ? X f  I  >  ? Z (  j  ■  •F.FTHa>«J-:*NNOn:- 

00300 

12 

F  ORMAT (3 *1 14*3F6.0)  ’ 

00310 

CAL..  DA?HLW(2.6  i 

0032) 

WRITE! 6? 15) 

00330 

WRITE (6? 16)  I J?X( J) ?Z( J) »BETA! J>  ? J-l ■  Hft : 

00340 

15 

£  Ok’tfAT  *  ■  1 6.X  * 4  { X  f  / )  COORD  [ rJA  1 L  AND  1 £  I C » l ON  COEF 

FOE  EACH  NOD " r 

00350 

1  / 1 X » 3  <  *  NOD  V  -KIC  •  /  ’ 

00  .-’60 

16 

FORMAT! 3( I4?2F8«1?F6*1) ) 

00370 

READ  !5  <  2)  !  J?  ( J  CON!  I  ?.J)  ?  I  - 1  •  7)  ?  1.  -  2  •  NELE ) 

00380 

CAtL  PASHLN(26) 

0038'  ■ 

WRITE! 6? 20) 

00400 

20 

FORMAT (/28>»* NODAL  CONNECTIVITY '/IX? 4! *ELEM  N1 

N2  f<3  * ;  / 

00410 

WRITE (6?21)  !l..?(  ICON(I?L  )?I-i*3'  ?L-I  .NELE' 

00420 

2! 

FORMAT! 4!4I5> ? 

00430 

CALL  DASHLN ■' Xe > 

00440 

NNBR=NNOD-NBK 

00 150 

DO  23  T-l  ,HW 

A/ 1  ,i  /-  /) 

l  N  ft  ( } )  -NWt’R^  t 

00470 

2"? 

CONT  INI  IF 

00180 

WRITE  1 6? 25 •  *BR 

•:w>49(-, 

25 

FORMAT! /12X?*  NUMBER  OF  NODES  ON  REFLEC-ION  DOMAIN?  NBR--’»I4/ 

00500 

3  27X  ?  *  THE  1  V !  iMNE  C  i  ’  V .  ;  V  ARE  1  •  > 

00519 

WRITE! 6 ?2)  (IMF-. I>? [-1.KPR, 

00520 

DO  26  T=1 r NRT 

00530 

INBT <  I 
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00540  26  CONTINUE 

00550  URITE(6»28)  NET 

00560  28  F0RHAT(//10X»*  NUMBER  OF  NODES  ON  TRANSMISSION  DOMAIN,  NBT=SI4/ 

00570  1  27X, ‘THEIR  CONNECTIVITY  ARE!') 

00530  WRITE(6»2)  <INBT<I),I=1,NBT> 

00590  WRITE<6»30)  NSGF 

00600  30  F0RMAT(//13Xf  ‘  NUMBER  OF  SEGMENTS  OF  FREE  SURFACE,  NS3F=M4) 

00610  DO  35  1=1 i NSGF 

00620  READ (5,2)  J J » < INF ( J » I > » J= 1 » J J > 

00630  NF(I)=JJ 

00640  WRITE<6,32>  I » I » NF < I ) 

00650  32  FORMAT (/14X»‘  NUMBER  OF  NODES  0N‘»I2,‘-TH  SEGMENT,  NF<SI2> 

00660  1  *)=‘,I4/27X, ‘THEIR  CONNECTIVITY  ARE  I  * ) 

00670  WRITE (6*2)  ( INF< J»I>, J=1,JJ> 

00680  35  CONTINUE 

00690  CALL  DA8HLN126) 

007 00  HR=-Z(NNBRH) 

00710  XR=X(NNBR+1) 

00720  HT=-Z(1) 

00730  XT=X ( 1 ) 

00740  WRITE 0 6 » 40 )  HR » XR * HT , XT 

00750  40  FORMAT </2X, ‘FOR  REFLECTION  DOMAIN! */13X» ‘  WATER  DEPTH.  HR=‘.F8.2» 

00760  1‘  HORIZONTAL  EXTENT,  XR=‘ ,F8.2»/2X, ‘FOR  TRANSMISSION  DOMAIN!*/ 

00770  2.1 3X, '  WATER  DE  PTH,  HT=*,F8-2,‘  HOR.ZONTAL  EXTENT,  XT=',F8.2/> 

00780  IF( IFORCE.EQ.O)  GOTO  50 

00790  CALL  DASHLNE26) 

00800  WRITE (6, 42)  NSGB 

00810  42  FORMAT <//7X»*  NUMBER  OF  BODY  FOR  FORCE  CALCULATION,  NS6B=  *,I4> 

00820  DO  45  1=1, NSGB 

00330  READ; 5, 2)  JJ, ( INBOD (J, I) ,J=1 ,JJ) 

00840  READ> 5,4 j  SCLE I ) ,XC( I s ,ZC(I> 

00850  NBODv  I  5=.JJ 

00860  WRITE (6, 44)  I,I,NBOD(I),SCL;i),XC(I),ZC(I) 

00870  44  FORMAT (/12X, ‘NUMBER  OF  NODES  0N‘,I2,*-TH  BODY,  NB0D<*,I2,‘ >=  ‘,14, 

00880  1  /30X ,  ‘LENGTH  SCALE,  SCL ( I )  =  * , F8 ♦ 2/25X , ‘ REFERENCE  CENTER,  <XC»ZC 

00890  2)=' ,2F8.2/6X, 'THEIR  CONNECTIVITY  ARE!*/ 

00 900  WRITE (6, 2)  <INB0D(J,I>,J=1,JJ> 

00910  45  CONTINUE 

00920  CALL  DASHLNE26) 

00930  50  CONTINUE 

00940  CALL  BAND ( I CON , NELE » NBR , NBT , NBD ) 

00950  400  CALL  DAT AI N < IGO , OMGSG , OMGA , 6 ) 

00960  IF (IGO. EG. 0!  GOT 0  404 

00970  C  CALL  DASHLN126) 

00980  CALL  EIGVAL ( WKR , WKHR , NKR , HR , OMGSG , P I ?  TOLR ) 

00990  CALL  EIGVAL ( WKT , WKHT , NKT , HT , OMGSG , PI , TOLR ) 

01000  C  CALL  DASHLN<26> 

01010  CALL  FRICTI<EPSN, XLAMD, BETA, NNODfOMGA, GAMA, VISCO.HR, WKR, PI) 

01020  CALL  LAPLAC ( 3YSK , ICON , NELE , NNOD , NBD , X , Z , XLAMD ) 

0 1030  CALL  SURFAC ( SYSK , NNOD , NBD , NSGF , NF , INF , NFMX , X , GhGSG ) 

01040  CALL  REFTRA (SY SDR, WKR, WKHF:, NKR, OMGSG »XR. 1.0) 

01050  CALL  REFTRA (SYSDT, WKT, WKHT, NKT, OMGSG, XT, -1.0) 

01060  CALL  HYBR I D < BYSKR , NNOD , WKR , WKHR , NKR , Z , XR , HR , I NBR , NBF: , 1 . 0 ) 

01070  CALL  HYBRIDISYSKT, NNOD, WKT, WKHT, NKT, Z, XT, HT.INBT, NBT, -1,0) 
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I 


01030 

01090 

01  too 
onto 
O’  120 
01130 
01 140 

j  1 160 
01170 
01 130 
Oil  ;'0 
01200 
0121 0 


CALL  CZERQ ( S  YSOR » NKR  > 

3YSQR <1)  =  <C*0»-1«0) MKR ( 1 )  <*AM  <  GHGSG » WKR ( 1 ) » WKHR ( 1 ) » 1 ) 

"ALL  CZERO ( SYSQT  >  NKT ) 

CALL  BOUNOA (SYSQ* NNOO * WKR (l)i UKHR' 1 ) , NBR*Z * XR » HR ) 

CALL  REDUCE (SYSK *SYSQ*SY5KR *SYSDR*SYSQft * NNGD*KBD *NKR*N8ri?2) 

CALL  REDUCE (SYSK , SYSQ * SYSKT * SYSHT*SYSG7 . NNGB , HBB * MKT * NET*  1 ) 

"ALL  CSIMQ (SYSK . SYSQ - i.'NOD.NBD ) 

CALL.  SOLVE (SYSQ * SYSKR , SYSDR * SYSOR » MNOD * NKR ,NBR*2> 
i ALL  SOLVE (SYSQ * S YSKT *SY3DT r SYSQT » NNOD . MKT *N8T?1 ) 

[F:,MM0D0T.l-.1>  GOTO  67 
WRITE (6? 60) 

60  FORMAT  (////20C  -’)/’  -  THE  SOLUTION  OF  THE  SYSTEM  - 

1  20 C  )/''27Y»  ’NODAL-POTENTIAL V12X* ’NODE  REAL-PART  IMAGE- PAr 

2T  ABS-VALUE  PHASE’/) 


01220 

DO  65  I  1-NNOD 

01230 

CALL  0U7CPX*  I.SYSQd)) 

01240 

65 

CONTINUE 

01250 

67 

IF (NKRUT.LI • i >  GOTO  77 

01260 

WRITE! f 70 i 

0  j  2 ' 0 

70 

FORMAT (,7/18X» ’COEFFICIENTS  FOR  REFLECTION  DOMAIN’) 

01280 

WRITE (6) 72  i 

01290 

72 

FORMAT (12X*  *COt'  REAL-PART  IMAGE-PART  ABS-VALUE 

01300 

DO  75  I  1. NT  ROT 

01310 

CALL  OUTCF'X ( I « SYSQR(  t ) > 

1 320 

75 

CONTINUE 

•.  S  yj 

/  / 

IF(NK'/V.  .1  T,i)  GOTO  S7 

01340 

WRITE <6- 80) 

01350 

SO 

FORMAT (///l 7X. ’COEFFICIENTS  FOR  TRANSMISSION  DOMAIN* 

01 360 

WRITE! 6 » "2 ' 

01370 

DO  85  I--1  »NKTO“ 

01380 

CALL  OUICP* ( I » SYSQT  *  I) ) 

01390 

35 

CONTINUE 

01400 

n  7 

IF* I FORCE . XL  • 0 >  CALL  FORCE ( SYSQ*  X • Z » NNOD  *NBuDMX*NSGB 

01410 

1  3CL *X  »ZC*tPSN*X L AMD * 0M3A 

01420 

CALL  DA6HLN(26) 

01430 

GOTO  400 

01440 

404 

CONTINUE 

01  45' 

CC 

CALL  CFflME(CB) 

01460  CC 
01470  CC 


PT-CB-CA 
WRITE* 6.500)  DT 


014E0  CCU-K)  FORMAT ( ,'9 CPU  TIME*  DT='*iPE12.4**  SECONDS  MO!’*’)) 

01490  STUI- 

01500  END 

01510  SUBROUTINE  DASHLN(N) 


01530 

01540 

01550 

01560 

01570 

01580 

01590 

01600 


DATA  DASH/3H-  -/ 

WRITE (6. 4)  (BASH*I=1*N) 
4  FORMAT* //IX. 42A3) 
RETURN 
END 

SUBROUTINE  ZERO(A*N) 


01620 
01630 
01640 
01650 
01660 
01670 
01680 
01690  C 
01700  C 
01710 
01720 
01730 
01740 
01750 
01760 
01770 
01780  C 
01790  C 
01800 
01810 
01820 
01630 
01840 
01850 
01860 
01870 
01880 
01890 
01900 
01910 
01920 
01930 
01940 
01950 
01960 
01970 
01980 
01990  C 
02000  C 
02010  C 
02020 
02030 
02040 
02050 
02060 
02070 
02080 
02090 
02100 
02110 
02120 
02130  C 
02140  C 
02150 


DIMENSION  A(N> 

DO  4  1=1, N 
A  < I ) =0 ♦ 0 
4  CONTINUE 
RETURN 
END 

SUBROUTINE  CZERO<A,N) 


COMPLEX  A(N) 

DO  4  1=1, N 
A(I)=<0. 0,0.0) 

4  CONTINUE 
RETURN 
END 

SUBROUTINE  BAND < ICON, NELE , NBR , NBT , NBD ) 


DIMENSION  IC0N(3,NELE> 

KMX=1 

DO  8  L=1,NELE 
IMX=ICONU,L> 

IMN=IC0N(1,L> 

DO  4  l<=2,3 

IRICON(K,L).LT.IMN>  IMN=ICON(K,L) 

IF ( ICON(K»L) iGT. IMX)  IMX=ICON(K,L> 

4  CONTINUE 
KK=IMX-IMN 

IRKRGT.KMX)  KMX=KK 
8  CONTINUE 
NBD=NMX+1 

NBD=AMAXO(NBD,NBR,NBT) 

WRITE<6, 10)  NBD 

10  F0RMAT(//27X, "BANDWITH,  NBD=" , I4//40( ")> 
RETURN 
END 

SUBROUT I NE  DAT AI N ( I GO , OMGSG , OMGA , G ) 


OMGSGRWAVE  FREQUENCY >**2/GRAVITY  CONSTANT, 


READ(5,4)  IGO,WAVT 
4  FORMATS  1 10, F10*l) 

IF(IGO.EQ.O)  RETURN 
WRITE (6, 8)  WAVT 

8  FORMAT ( //20 ( 1  — '  )/3X,  'WAVE  PERIOD,  WAVT=" ,F8.2,  *  SECONDS"  ,/20*. 

1  )//> 

0MGA=6,2831853/WAUT 

OMGSG=OMGA*OMGA/G 

RETURN 

END 

SUBROUTINE  E 1 6UAL ( WK , WKH , NK , H , OMGSG , P I , T OLR ) 


DIMENSION  WN(NK')  ,WKH(NK) 
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02160  C 
02170  C 
02180  C 
02190 
02200 
02210 
02220 
02230 
02240 
02250 
02260 
02270 
02280 
02290 
02300 
02310 
02320 
02330 
02340 
02350 
02360 
02370  C 
02380  C 
02390  C 
02400  C 
02410  C 
02420  C 
02430  C 
02440 
02450 
02460 
02470  C 
02480  C 
02490 
02500 
02510 
02520 
02530 
02540 
02550 
02560 
02570 
025B0 
02590 
02600 
02610  C 
02620  C 
02630 
02640 
02650 
02660 
02670 
02680 
02690 


WRITE (6. 4) 

4  r0RHAT(//45<’-,)/4<,-*)f '  EIGEN  VALUES  OF  DISPERSION  EQUATION  ** 
14( )/45( '-’)//19X» "FIRST  NUMBER  IS  REAL*  OTHERS  ARE  IMAGINARY’/) 
C=OMGSG*H 
XJ=C 
10  XI=XJ 

XJ=C/TANH(XI> 

IF<ABS< XI-XJ) .GT.TGLRfcXJ)  GOTO  10 
WKH< 1 )=XJ 
WK<1)=WKH<1)/H 
IF<NK<LE<1 )  GOTO  25 
DO  20  1=2 *NK 
XJ=(I-1)*PI 
DX=X.J 
15  Xl^XJ 

X  J=AT AN ( -C/XI ) +DX 
IF<ABS<XI-XJ).GT,TOLR*XJ)  GOTO  15 
WK'H<  I  )=XJ 
UK<  J )-UKH<I)/H 
20  CONTINUE 
25  CONTINUE 
WRITE (6 *35) 

35  FORMAT <1 IX >’ WAVE  NUMBER*  UK*  ARE 5’) 

WRITE<6*37)  <WK<I>*I=1*NK) 

37  FORMAT <8F1 0.4) 

WRITER* 40) 

40  FORMAT ( '  WAVE  NUMBER  *  DEPTH*  WKH*  ARE i ' ) 

WRITE<6*37)  (WKH(l>.I=.t*NK) 

RETURN 

END 

SUBROUTINE  FRICTI <EPSN*XLAMD*BETA*NNOD*OMGA*GAMA*VISCOjH*WK*FI ) 


COMPLEX  EPSN < NNOD ) *  XLAMIK NNOD ) * TM* WI * EXP3 
DIMENSION  BETA (NNOD) 

WI-( 0.0*1 .0>/0MGA 

EXF'G=COS  ( GAMA  )  +  <  0 , 0 . 1 . 0  >  #SI N  <  GAMA ) 

TM=EXF’G*SQRT(  VISCOIOMGA) 

DO  20  1=1. NNOD 

EPSN ( I ) =TM*BET A(I) 

XLAMD<I)=1.0/(1.0+WI*EPSN{I)) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  EIGVAL1(WKH*H*0MG3G*T0LR) 


C=OMGSG*H 
XJ=C 
10  XI=XJ 

XJ=C/TANH(XI) 

IF (ABS( XI-XJ ) <GT*  TOLRIXJ)  GOTO  10 

WKH=XJ 

RETURN 
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02700 
02710 
02720  C 
02730  C 
02740  C 
02750 
02760 
02770 
02780 
02790 
02800 
02810 
02820 
02830 
02840 
02850 
02860 
02870 
02B80 
02890 
02900 
02910 
02920  C 
02930 
02940 
02950  C 
02960  C 
02970 
02980 
02990 
03000 
03010  C 
03020  C 
03030  C 
03040 
03050 
03060 
03070 
03080 
03090 
Oil  00 
03110 
03120 
03130 
03140 
03150 
03160 
03170 
03180 
03190 
03200 
03210 
03220 
03230 


END 

SUBROUTINE  ELMKIMEL.X, Y.ELK) 


GENERATION  OF  TIANGULAR  ELEMENT  MATRIX.  ELK. 


DIMENSION  X(3),Y<3>»ELK{3»3) 

CALL  ZERO (ELK. 9) 

B1=Y(2)-Y(3) 

B2=Y(3)-Y< 1) 

B3=Y(1)-Y(2) 

C1=X(3)-X(2> 

C2-X(l)-X(3) 

C3=X(2)-X< 1 ) 

AREA=0.5*(B1*C2-B2*C1) 

IF ( AREA. GT. 0.0)  GOTO  20 
WRITE(6.10)  MEL. AREA 

10  FORMAT C  #****» 15, *-TH  ELEMENT,  AREA-' »F10.5) 

STOP 

20  A4=4. KAREA 

ELK(1,1)=(BI*B1+C1*C1)/A4 
£LK<  1  »2)=-(Bl#B2+Cl*C2)/A4 
ELK ( 1 ,3)- (B1*B3+C1*C3)/A4 
ELK(2,1)=ELK(1»2) 

ELK ( 2 , 2 )  =  < B2*B2+C2*C2 ) /A4 
ELK ( 2 , 3 )  =  ( B2*B3  f C2*C3 ) /A4 
£LK(3»1)=ELK(1»3) 

ELK(3,2)=ELK(2»3) 

FLM3,3)  =  (B3*B3+C3*C3)/A4 

RETURN 

END 

SUBROUT  I NE  L  AF'L  AC  <  S  YSK ,  I  CON ,  NELE ,  NNOD ,  NBD » X ,  Y ,  XL  AMD  » 


ASSEMBLE  ELK  INTO  SYSK  FOR  LAPLACIAN. 


COMPLEX  S YSK ( NNOD , NBD ) , XL AMD ( NNOD  > , TM 

DIMENSION  I CON ( 3 , NELE ) » X ! NNOD ) » Y ( NNOD ) » I E ( 3 ) , XE ( 3 ) , YE ( 3 > , ELK ( 3 » 3  > 

NN=NNOD*NBD 

CALL  CZERO(SYSK.NN) 

DO  40  L=1 .NELE 
TM=<0. 0.0,0) 

DO  10  J-l .3 
IE( J)=ICON( J.L) 

XE( J)=X(IE(J> ) 

YE( J)=Y( IE( J> ) 

TM=TM+XLAMD< IE( J) ) 

10  CONTINUE 
TM=TM/3. 

CALL  ELMK(L»XE» YE.ELK) 

DO  30  1=1,3 
DO  20  J=I,3 
JI=IE(J)-IE(I) 

IF(JI.GE.O)  GOTO  15 
JI1=-JI+1 

SYSK(IE(J),JI1)=SYSK(IE(J)»JI1)+TM*ELK(I,J) 
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03240 

03250 

03260 

03270 

03280 

03290 

03300 

03310 

03320 

03330  C 

03340  C 

03350 

03360 

03370 

03380 

03390 

03400 

03410 

03420 

03430 

03440 

03450 

03460 

03470 

03480 

03490 

03500 

03510 

03520 

03530 

03540 

03550 

03560 

03570 

03580  C 

03590  C 

03600  C 

03610  C 

03620 

03630 

03640 

03650 

03660 

03670 

03630 

03690 

03700 

037.10 

03720 

03730 

03740 

03750 

03760  C 

03770  C 


GOTO  20 
15  JI1=JI+1 

SYSKUEU)*  JIl)=SYSK(IE(I)f  JIi)+TM*ELK(I*J> 

20  CONTINUE 
30  CONTINUE 
40  CONTINUE 
RETURN 
END 

SUBROUT I NE  SURFAC (SYSK* NNOD , NBD * NSGF * NF * INF * NFMX . X * OMGSG ) 


COMPLEX  SYSK (NNOD* NBD) 

DIMENS  I ON  X ( NNOD  > . NF ( NSGF ) *  INF < NFMX * NSGF ) 
EX=0MGSG/3. 

DO  30  1=1* NSGF 
J J=NF ( I ) 

DO  20  J=2» JJ 
J1=INF(J-1*I> 

J2=INF< J*I> 

EXL=EX*(X< J2)-X( Jl)) 

S YSK ( J 1 > 1 ) =S YSK < J 1 * 1 ) -EXL 
S  YSK ( J2 . 1 ) =SYSK ( J2  *  1 ) -EXL 
J12-J1-J2 

IF( J12.GE.0)  GOTO  10 
JA--J12+1 

S YSK ( Jl * JA ) =SYSK ( Jl* JAY -EXL  72 . 

GOTO  20 
10  JA— J 12+1 

S YSK ( J2  j  JA ) -SYSK ( J2*JA) -EXL/2 . 

20  CONTINUE 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  REFTRA(SYSD*WK*WKH*NK* OMGSG? XX* SGN ) 


ASSEMBLE  ELEMENT  OF  COEFFICIENT  TYPE. 

USE  SGN-1 . 0  FOR  REFLECTION.  SGN=-1.0  FOR  TRANSIMISSION. 


COMPLEX  SYSD(NK) »CX 
DIMENSION  UK(NK) .WKH(NK) 

CALL  CZERO(SYSD.NK) 

DO  20  1=1* NK 
XK=SGN*2.*UK(I)*XX 
IF(I.GT.l)  GOTO  10 
CX=(0.0*I.0)*COS(XK)-SIN(XK) 

SYSD ( I ) =CX*UK ( I ) *AM  ' OMGSG *  UK  < I ) *  WKH ( I ) *  I ) 

GOTO  20 

10  SYSD(I)=-EXP(-XK)*WK(I)*AM(OMGSGfWK(I),WKH(I)»I) 
20  CONTINUE 
RETURN 
END 

FUNCTION  AM ( OMGSG  >  WK  *  WKH  *  I NDX ) 


INTEGRATION  OF  (C0SHK(Z+H)/C0SHKH>**2  FROM  -H  TO  0. 
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03780  C 

03790  C 

03600  C 

03810  C 

03820  C 

03830  C 

- 

03840 

03850 

03860 

03870 

4 

03880 

03890 

03900 

03910  C 

-  ■ 

03920  C 

03930  C 

-  • 

03940 

03950 

03960 

03970 

03930 

03990 

04000 

04010 

04020 

04030 

04040 

04050 

10 

04060 

14 

04070 

04C80 

04090 

04100 

04110 

04120 

04130 

24 

04140 

28 

04150 

04160 

04170 

30 

04180 

04190 

04200 

40 

04210 

04220 

04230 

04240  C 

-  - 

04250  C 

-  - 

04260 

04270 

04280 

04290 

04300 

04310 

OMGSG=OM*OM/G.  WHERE  OM=WAVE  CIRCULAR  FREQUENCY  AND  6=GRAVITATI0N 
CONSTANT. 

WK:=WAVE  NUMBER 

WKH=WK*H.  WHERE  H=WATER  DEPTH 

INDX=lf  PROPAGATING  MODE.  OTHERWISE  EVANESCENT  MODES. 


IF<INDX.EG.1>  GOTO  4 

AM=0 . 5*  < WKH/COS ( WKH ) **2-0MGSG/WK ) /WK 

RETURN 

AM=0 . 5* < WNH/COSH ( WKH  5  **2+QM6SG/WK ) /WK 

RETURN 

END 

SUBROUT I NE  HYBRID ( SYSKRT . NNOD .WK.WKH.NK.Z.XX.H. INB.NB.SGN) 


USE  SGN=1 .0  FOR  REFLECTION.  SGN=-1.0  FOR  TRANSIMISSION. 


COMPLEX  SYSKRT ( NK . NB ) . CX . TM 
DIMENSION  Z(NNOD) »WK(NK) »WKH(NK) »INB(NB> 

NN=NK*NB 

CALL  CZERO< SYSKRT. NN) 

DO  40  1=1 .NK 
EX=SGN#WK(I)*XX 
WKZ2=0.0 
CH2=1 ,0 

IF(I.EQ.l)  GOTO  10 
CX=EXP(-EX)/COS<WKH(I)) 

GOTO  14 

CX=<  (0.0.1 .0)*C0S(EX)-SIN(EX) )/COSH(WKH(I ) ) 

DO  30  J=2.NB 
WKZ1=WKZ2 

WKZ2=WK(I)*(Z(INB(J))+H) 

CH1=CH2 

IF(I.EQ.l)  GOTO  24 
CH2=C0S(WKZ2) 

GOTO  28 
CH2=C0SH(WKZ2) 

TM=CX*(CH2-CH1 )/<WKZ2-WKZl ) 

SYSKRT ( I . J-l ) =S YSKRT  < I . J-l ) -TM 
SYSKRT (I. J) =S YSKRT ( I . J ) +TM 
CONTINUE 

IF  < I . EQ . 1 )  SYSKRT ( I . NB >  =S YSKRT ( I .NB  > -CX*SINH( WKH ( I > ) 

IF ( I . NE . 1 )  SYSKRT ( I , NB ) =S YSKRT ( I . NB ) +CX*SIN ( WKH ( I > ) 

CONTINUE 

RETURN 

END 

SUBROUT  I  NE  BOUNDA  ( S  YSO .  NNOD .  WK  1 .  WKH1 .  NB » Z » XX .  H ) 


COMPLEX  SYSQ(NNOD) .CX.TM 
DIMENSION  Z(NNOD) 

CALL  CZERO(SYSQ.NNOD) 

NN=NNOD-NB 

XK=-WKltXX 

CX=((0.0»1.0)*C0S(XK)-SIN(XK))/C0SH(WKH1) 
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04320 

04330 

04340 

04350 

04360 

04370 

04380 

04390 

04400 

04410 

04420 

04430 

04440 

04450 

04460 

04470 

04480 

04490  C 

04500  C 

04510  C 

04520  C 

04530 

04540 

04550 

04560 

04570 

04580 

04590 

04600 

04610 

04620 

04630 

04640 

04650 

04660 

04670 

04680 

04690 

04700 

04710 

04720 

04730 

04740 

04750  C 

04760  C 

04770 

04780 

04790 

04800 

04810 

04820 

04830 

04840 

04850 


WKZ2=0,0 
CH2=1.0 
DO  10  J=2»NB 
NJ=NN+J 
WKZ1=WKZ2 

UKZ2=WK1*(Z(NJ)+H) 

CH1=CH2 

CH2=C0SH(WKZ2) 

TM=CX*  < CH2-CH1 ) / < WKZ2-WKZ1 ) 

SYSQ(NJ-1)=SYSG(NJ-1)-TM 
SYSQ(NJ)=SYSO(NJ)+TM 
10  CONTINUE 

SYSQ ( NNOD ) =SYSQ <  NNOD ) -CX*SI NH ( WKH1 ) 

RETURN 

END 

SUBROUTINE  REDUCE < SYSK » SYSQ » SYSKRT , SYSDRT * S YSQRT r NNOD » NBD » NK * NB > 
1  INDX) 


INDX=1»  MERGE  FROM  THE  TOP  OF  THE  MATRIX i  OTHERWISE  MERGE  FROM  THE 
BOTTOM. 


COMPLEX  SYSK ( NNOD » NED  ), SYSQ ( NNOD ) > SYSKRT <  NK » NB ) > SYSDRT <  NK )  > 
1  SYSGRT (NK) •TMjTMOjTMI 

NN=NNOD-NB 
DO  30  I-1..NK 
TMO -SYSGRT ( I ) /SYSDRT ( I ) 

DO  20  J=1»NB 
J2=-J 

IF ( INDX.NE . 1 )  JJ=NN+J 
SYSG(JJ)=SYSG(JJ)-TMO*SYSKRHI»J) 

TM1=SYSKRT(I> J)/SYSDRT(I) 

DO  10  K=1 ? J 

KK"K 

IFdNDX.NE.l  >  KK=NN4K 
JK=J-N+1 

TM=TMl*SYSKRT(IfK) 

SYSK ( KK  t  JK ) =SYSK { KK » JK ) -TM 
10  CONTINUE 
20  CONTINUE 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  CSIMQ<A»B*NEGT»NBD) 


COMPLEX  ArB 

DIMENSION  A(NEQT»NBD) »B(NEQT) 

T0LR=1 .0E-30 

NNB=NEGT-NB0+1 

NB=NBD-1 

DO  24  I=1»NEQT 

IF(CABS( A( I>1) ) .LE.TQLR)  GOTO  40 
A(If 1)=1.0/A(I>1) 


BIO 


04860  IF(NNB.GE.I)  GOTO  10 

04870  NB=NEQT-I 

04880  IF(NB.LE.O)  GOTO  24 

04890  10  Nl=l 

04900  BO  14  J=1fNB 

04910  Nl-Nl+1 

04920  14  A  < I » N1 ) =A  < I » N1 ) (1*10 

04930  BO  20  J=1 fNB 

04940  IJ=I+J 

04950  JJ=NB+1-J 

04960  J1=J+1 

04970  Baj)=B(IJ)-B(I)*Aa»Jl)/Aa»l) 

04980  DO  16  K=1.JJ 

04990  KJ=K+J 

05000  16  A<IJfK)=A(IJ»K)-A(I»J1)#A{I»KJ)/A(If1) 

05010  20  CONTINUE 

05020  24  CONTINUE 

05030  J=NEQT 

05040  Nl-NEQT-1. 

05050  DO  34  K=1fN1 

05060  NB-NBD-1 

05070  J=J-1 

05080  IF(NNB.GE.J)  GOTO  28 

05090  NB=NEQT-J 

05100  28  DG  34  1=1. NB 

05110  IJ=J+I 

05120  JJ=I+1 

05130  B<J)=B(J)-B<IJ)*A(JfJJ) 

05140  34  CONTINUE 

05150  RETURN 

05160  40  WRITE<6.41)  A<If1>.I 

05170  41  FORMAT ( '  *  ELEMENT  OF  DIAGONALS f2E15.6f *  AT*fI4f*-TH  ROW  *  ■) 

05180  STOP 

05190  END 

05200  SUBROUTINE  SOLVE ( SYSQ  f  S YSKRT f  S YSDRT  f SYSQRT f  NNOD f NK f NB f INDX ) 

05210  C - - - 

05220  C  INDX=1f  MERGE  FROM  THE  TOP.  OTHERWISE  FROM  THE  BOTTOM, 

05230  C  - - - 

05240  COMPLEX  SYSQ ( NNOD ) . SYSKRT ( NK f  NB ) f S YSDRT ( NK ) f  SYSQRT < NK ) . TM 

05250  NN=NNOD-NB 

05260  DO  20  1=1 fNK 

05270  DO  10  J=1.NB 

05280  JJ=J 

05290  IF(INDX.NE.l)  JJ=NN+J 

05300  TM=SYSQ(JJ)*SYSKRT(I.J) 

05310  SYSQRT ( I)=SYSQRT(I )-TM 

05320  10  CONTINUE 

05330  SYSQRT ( I ) =S YSQRT < I ) /SYSDRT ( I ) 

05340  20  CONTINUE 

05350  RETURN 

05360  END 

05370  SUBROUTINE  OUTCPX(IfCA) 

05380  C - - - - 

05390  C . . . . . 
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05400 

05410 

05420 

05430 

05440 

05450 

05460 

05470 

05480 

05490 

05500 

05510  C 

05520  C 

05530 

05540 

05550 

05560 

05570 

05580 

05590  C 

05600  C 

05610  C 

05620  C 

05630  C 

05640 

05650 

05660 

05670 

05680 

05690 

05700 

05710 

05720 

05730 

05740  C 

05750 

05760 

05770  C 

05780 

05790 

05800 

05810 

05820 

05830 

05840 

05850 

05860 

05e70 

05880 

05890  C 

05900 

05910 

05920  C 

05930 


COMPLEX  CA 
AS -REAL (CA) 

AI -AIMAG(CA) 

AT-ATAN2<AI»AR) 

AS--SGRT<AR*AR+AI*AI) 

WRITE(6r4)  I »ARf AI • ASt AT 
4  FORMAT (8Xf I7»1X»3F12, 4fF7.3) 

RE  TURN 
ENK 

SUBROUTINE  FORCE  <  SYS9 » X i Z  > NNQD , NBODMX » NS8B  >  NBOB » INBOB  *  3CL » XC » ZC » 
1  EPSN  > XLAMB i OMGA ) 


COMPLEX  SYSQ(NNOD) f TM»P1 »P2*FX»FZ» VMY*EF'SNf NNOB) »XLAMJ;<  fJWQD) .Tit 
1  T2»TXf TZrWI 

DIMENSION  X(NNOD) >Z(NNOD) fNBOD(NSGB) > INBOD ( NBODHX»NSGB )» SCL (NS6B ) ? 
1  XL'(NSGB)  f  ZC(NSGB) 

WI=(0.0>1.0)/0MGA 
DO  30  1=1 eNSGB 
WRITE < 6 » 5 )  I 

5  FORMAT ( ///l 9X»’ HYDRODYNAMIC  FORCES  F0RM2»*-TH  BODY') 

WRITE' 6r 10) 

10  FORMAT (//31X> ’PRESSURE* /12X> ‘NODE  REAL-PART  IMAGE-PART  ABS-VA 
1LUE  PHASE*/) 

JJ=NBOD(I ) 

FX=(0,0r0.0) 

FZ=(0. 0.0.0) 

VMY=<0. 0,0.0) 

J1=INB0D(JJ»I) 

J2=INB0D( 1 > I ) 

X1=X(J1)-XC<I) 

Z1=Z<J1)-ZC(I) 

X2 -X  <  J2 ) — XC ( I > 

Z2-Z(J2)-ZCA1) 

AREAB=-0.5*(X2-X1)*(Z1+Z7) 

T2=-WI *EPSN ( J2 ) *XL AMD ( J2 ) 

P2=SYSQ(J2) 

CALL  OUTCPX( J2»P2) 

DO  20  J=2* JJ 

J1=J2 

X1=X2 

Z1=Z2 

F'1=P2 

T1=T2 

J2=INBGB< J»I) 

X2=X( J2)-XC< I ) 

Z2=Z( J2)-ZC(I) 

P2=SYSQ(J2) 

T2=-WI*EPSN<J2)*XLAMB<J2) 

CALL  0UTCPX( J2»P25 

DX=X2-X1 

BZ=Z2-Z1 

AREAB=AREAB-0 . 5*DX* (Z1FZ2) 

TM=0.5*(P1+P2) 
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05940 
05950 
05960 
05970 
05980 
05990 
06000 
06010 
06020 
06030 
06040 
06050 
06060 
06070 
06080 
06090 
06100 
06110 
06120  C 
06130  C 
06140  C 
06150 
06160 
06170 
0&180 
06190 
06200 
06210 
06220 
06230 
06240 
06250 
06260 


FX=FX~TM#DZ 

FZ=FZ4TM*bX 

TH=0,5t(TlU2)*(P2~Fl) 

TX=TM*BX 

TZ=TM*DZ 

FX=FX+TX 

FZ=FZ+TZ 

X11=X1*X1 

X12=X1*X2 

X22=X2*X2 

Zll^ZlIZl 

Z12=Z1*Z2 

Z22=Z2*Z2 

TM--P1*(Z22+Z12-2.0*Z11)+P2*(2.0«22-Z12-Z11) 

TM=(TM+Pl*<X22+X12--2.0*Xll)+P2*<2.0*X22-Xi2-Xll>>/6.0 

WfY=VMY+Trt 

0MY-VMY4 < -TX*  <  Z24Z1 ) 4TZ*<  X24X1 >1/2,0 

20  CONTINUE 
FX-FX/AREAB 
FZ=F2/AREAB 
UMY~WY/AREAB/SCL(I ) 

FX--KX/SCl.’.I) 

FZ=FZ/SCl. ( I) 

OMY=OrtY/SCL(l>/SCL(I> 

WR1TE<6»24) 

24  FORMAT (//19X> 'WAVE  FORCES  1  INDXd=FX,2=FZ>3=MY  */12X.  *  INBX  REAL 
1ART  IMAGE-FART  ABS-VALUE  PHASE*/) 

CALL  0UTCF'X(1,FX) 

CALL  0UTCPX(2fFZ) 

CALL  0UTCF'X(3>VMY) 

30  CONTINUE 
RETURN 
END 
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APPENDIX  C:  NOTATION 

aQ  Incident  wave  amplitude 
A  Water  region 

B  Solid  body 

Wave  force  coefficient 
F  Total  force 

g  Gravitational  acceleration 

h  Water  depth 

i  Y-f 

j  subscript  J  =  1 , 2  represents  the  x  and  z  components 
+ 

km  Wave  number  of  the  evanescent  modes  m  ^  1 
+ 

k"  Wave  number  of  the  propagating  mode 

1  Length 

L  Length  scale  for  the  force  normalization 

m  Moment  coefficients 

M  Moment 

n  Unit  normal  vector  outward  from  the  water  region 

p  Pressure 

Pd  Hydrodynamic  pressure 

?  Distance  vector 

R  Coefficients  in  Equation  14 

t  Temporal  coordinate 

u  Flow  velocity 

U  Spatial  part  of  the  flow  velocity 

x  Horizontal  direction 

z  Vertical  direction 

a  Absorption  coefficient 

B  Friction  coefficient 

y  Phase  difference 

e  Friction  coefficient  of  the  solid  or  structural  boundaries 

C  Spatial  part  of  the  free  surface  displacement 

n  Free  surface  displacement 


Cl 


X  Associated  with  bottom  friction  and  wave  frequency  given  in  Equation  6 
v  Water  viscosity 

n  Functional 

p  Water  density 

t  Friction 

t  Velocity  potential  function 

•  Spatial  part  of  the  velocity  potential  function 

u>  Radian  wave  frequency 

3  Boundary  curve:  such  as  3A  ,  3B  ,  3F,...etc. 

V  Two-dimensional  gradient  operator 

o  Subscript  o  indicates  incident  wave 

♦  Superscript  +  indicates  the  reflection  water  region 

Superscript  -  indicates  the  transmission  water  region 


C2 


r 


